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A METHOD OF ESTIMATION 

The present invention relates to the field of non-invasive 
determination of bone quality of vertebrates. This determina- 
tion may be used in the diagnosis of, e.g., osteoporosis and 
5 other bone diseases which cause bone fragility and, thus, 
increase the risk of bone fracture. 

Accelerated bone loss leading to osteoporosis is a common 
phenomenon in women after menopause . Women often ignorantly 
suffer from accelerated bone loss, and the reduced strength 

10 of the bones is not discovered until a bone is broken or a 
vertebra collapses due to a load which a healthy bone or 
vertebra should be able to withstand. Thus, a large part of 
the osteoporotic patients are not aware of the reduced 
strength of their bones until a fracture reveals the disease 

15 and the extent thereof. 

Today, osteoporosis affects more than one third of elderly 
women in the industrialized part of the world. The prevalence 
of this disease is still increasing, partly caused by the 
increase in the proportion of elderly people, partly for 

20 unexplained reasons. In the lesser developed parts of the 
world, like South-east Asia and South America, it is pre- 
dicted that osteoporosis will become an enormous socio-econo- 
mic burden within the next 20-30 years. However, if indivi- 
duals at risk of developing osteoporosis can be identified, 

25 preventive measures can be applied. This requires a reliable, 
cheap and safe method for identification of those at risk. 

Since the risk of osteoporosis is closely related to the bone 
resistance to fracture, bone strength or bone quality mea- 
surements are likely to supply the essential information. 

30 Since the means of preventing osteoporosis are much more 

efficient than those of treating osteoporosis, identification 
of individuals at risk of developing this disease is crucial. 
Only by early prevention of osteoporosis, the individual as 
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well as the socio-economic consequences of the disease may be 
minimized. 



The method according to the invention provides a method for 
detection of reductions in bone quality, such as architecture 
5 and strength, at an early stage. Thus, the method of the 

invention provides a tool of screening potential patients for 
bone diseases and, thus, provides the possibility of early 
detection of bone diseases at a stage in which no symptoms of 
disease are noted by the patient. 



10 A typical measure of bone strength has been the Bone Mineral 
Density (BMD) of the bone. BMD measurements are typically 
obtained by X-ray of a bone together with a standard wedge. 
Having determined which part of the wedge attenuates the X- 
ray beam to the same degree as the illuminated bone, a 

15 measure of BMD may be obtained. 



A more sophisticated method of determining BMD is by X-ray 
absorptiometry of two different wavelengths. Using two wave- 
lengths enables the method to compensate for the effects of 
soft tissue etc. around the bone and, thus, to obtain a more 
20 exact determination of the X-ray attenuation of the bone. 

However, as will be clear from the following, the Bone Mine- 
ral Density of a bone is not necessarily connected to the 
actual strength of the bone. The reason for this is to be 
found in the structure of a bone of a vertebrate. 



25 A bone of a vertebrate consists of a cortical outer layer and 
a cancellous inner structure. Omission of the cancellous 
inner structure of a bone would result in a quite fragile 
bone. The cancellous inner structure of the bone consists of 
so-called trabeculae. Thicker vertical trabeculae are 

30 positioned in the bone in the direction of the main load 
(main compression or pull in the bone} , and thinner, 
horizontal trabeculae interconnect the vertical trabeculae. 
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Thus, the main density of a bone is constituted by the corti- 
cal layer and the vertical trabeculae. This is quite natural, 
as the bone reserves the largest part of the material to 
withstand the most common loads* A BMD measurement 
5 determining the density of a bone will therefore primarily 
estimate the amount of bone in the cortical layer and the 
vertical trabeculae and, thus, only to some degree the 
ability of the bone to withstand the loads to which the bone 
is adapted. 

10 However, when loads are applied not in the direction of the 
vertical trabeculae, the structure of the bone would be 
fragile without the thinner, horizontal trabeculae which 
interconnect the vertical trabeculae. The horizontal 
trabeculae define the ability of the bone to withstand loads 

15 not in the direction of the vertical trabeculae. Thus, the 
strength of the bone in directions not in the direction of 
the vertical trabeculae is defined by a very small part of 
the bone mass. 

Different bones vary in the relation between cortical and 
20 trabecular bone structures. Furthermore,, the architecture of 
the trabeculae vary according to the potential loading of 
this part of the bone. The number and diameter of 
longitudinal trabeculae and the presence of transverse 
trabeculae are also of fundamental importance for the 
25 fragility of the bone. These differences are of major 

importance for the strength of the bone, but they can only to 
a minor extent be detected by BMD measurements. 

Since trabecular bone has a greater surface area/weight 
ratio, it is to a higher extend exposed to accelerated bone 

30 loss than the cortical bone. Therefore, osteoporotic frac- 
tures are mainly seen in bones with predominantly trabecular 
structure, such as the distal forearm, proximal humerus, 
thoracic and lumbar spine as well as the femoral neck. 
Changes in BMD related to loss of trabecular bone severely 

35 underestimates the loss of resistance to fracture. As a 
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consequence, BMD measurements of patients with osteoporotic 
fractures reveal a major overlap with BMD measurements of 
normal subjects of the same cohort. 

The above information explains why persons with low energy 
5 fractures may have normal BMD values when compared to age- 
matched normal controls. Thus, one may wish to distinguish 
between bone quantity and bone quality, where the bone qual- 
ity is more closely related to the mechanical or bio-mechan- 
ical strength of the bone. 

10 A measure of the overall bone strength may naturally be 
obtained from a bone specimen taken from the potential 
patient and subjected to mechanical testing. However, this 
requires bone biopsy, which is painful and implies a minor 
risk of complications. Thus, in order to have a comfortable, 

15 cheap, fast and safe screening of the very large group of 
potential patients (most women after menopause) , the estima- 
tion of the bone quality should be performed non-invasively 
in a safe and fast manner. 

The present invention offers a non- invasive method which 
20 provides measures of the overall strength of bones and which 
is safe, comfortable for the potential patient, fast and 
cheap . 

In "A New Method for Automatic Recognition of the Radio- 
graphic Trabecular Pattern" Wil G.M. Geraets et al . , Journal 

25 of Bone and Mineral Research, Vol. 5. No. 3, 1990, pp 227- 
233, a method of recognizing the trabecular pattern from X- 
ray pictures is disclosed. According to this method, two 
types of noise are removed from an x-ray picture which has 
been scanned into an image memory: high-frequency noise is 

30 removed by using a median filter, and low- frequency noise is 
removed by using a local averaging operation on the image. 

After noise reduction a segmentation of the image data is 
performed whereby the image is binarized and subsequently 
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eroded (by using a median axis transform) so as to retain 
only central lines having a thickness of 1 pixel . A total of 
7 features are derived from the segmented and the eroded 
images. These features are subsequently correlated to the 
5 Bone Mineral Density measurements of the bones. 

However, due to the almost insignificant correlation between 
the image features obtained and the BMD measurements in the 
above-mentioned reference, it is evident that a higher cor- 
relation and, thus, a more precise estimate of the bone 
10 quality is required. 

In "computerized Radiographic Analysis of Osteoporosis: 
Preliminary Evaluation", Plilip Caligiuri et al, Radiology 
1993, 186; 417-474, a method in which X-ray images of the 
lumbar spine are scanned analyzed. In this reference, the 

15 power spectrum was obtained and the Root Mean Square (RMS) 
and the First Moment (FMO) thereof were compared to typical 
BMD measurements performed on the bones. It was concluded 
that both BMD did not correlate well with FMO and RMS and 
that this may be due to FMO and RMS containing additional 

20 information. In this reference, only the power spectrum and 
features derivable therefrom were investigated. 

The present invention provides a method and a system for 
providing an easily performed, but highly reliable estimation 
of the bone quality. 

25 In a first aspect, the present invention relates to a method 
for estimating the bone quality of a vertebrate, on the basis 
of two-dimensional image data comprising information relating 
to the trabecular structure of at least a part of a bone of 
the vertebrate, the image data being data obtained by expos - 

30 ing at least the part of the bone to electromagnetic radi- 
ation, the method comprising subjecting the image data to a 
statistical analysis comprising 
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a background correction procedure in which low frequent 
intensity variations not related to the trabecular struc- 
ture of the bone is reduced relative to image data 
related to the trabecular structure of the part of the 
bone, 

an image manipulation and feature extraction procedure 
wherein at least the local image intensity information as 
well as variation in the local intensity are utilized to 
extract information related to the trabecular structure 
of the part of the bone, the image manipulation and 
feature extraction procedure comprising subjecting the 
image data to at least one of the following procedures: 

(a) obtaining an estimate of the parametric estimate 
of the power spectrum of the image data and 
extracting features relating to the energy 
distribution of the parametric estimate, 

(b) obtaining, on the basis of image data on which a 
Fourier method has been used to emphasize the 
information in the image data relating to the 
trabecular structure, an estimate of a grey-level 
co-occurrence matrix and extracting at least one 
feature on the basis of the estimated co-occur- 
rence matrix, 

(c) obtaining an estimate of the projected trabecular 
pattern of the image data by using a Fourier 
method to emphasize the information in the image 
data relating to the trabecular structure and 
subjecting the manipulated image data to a mor- 
phological operation, and extracting features 
relating to the trabecular structure from the 
estimated projected trabecular pattern, 

(d) obtaining, on the basis of a frequency analysis 
of the image data, features relating to the 
periodicity of the trabecular structure of the 
part of the bone, 

and an estimation procedure in which the bone quality of 
the vertebrate is estimated on the basis of the derived 
features and optionally other features related to the 
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bone or the vertebrate and a predetermined relationship 
between the values of such features and reference bone 
quality parameters. 

As indicated above, in the present context "bone quality" is 
5 not equalled to "bone quantity", such as Bone Mineral Den- 
sity, as even a small loss of bone quantity may lead to a 
large loss of bone quality if the bone loss has taken place 
at critical positions in the trabecular structure. Thus, in 
the present context "bone quality" is a measure closely 
10 related to the risk of fracture of the bone, as it has been 
demonstrated in stress/strain evaluations of bone biopsies in 
which the microscopic structure is also evaluated (See, e.g., 
Lis Mosekilde relating to this issue) . 

Even though the bones and, thus, the trabecular structures 
15 are inherently three-dimensional, a projection of this struc- 
ture into two dimensions, such as a radiographic image, 
yielding a so-called projected trabecular pattern, conveys 
sufficient information about the trabecular structure and, 
thus, the bone quality to give a satisfactory estimate of the 
20 bone quality. 

Information relating to the trabecular structure may be local 
variations in, e.g., grey level information or any other 
information from which features relating to the trabecular 
structure may be derived. 

25 Low frequency intensity variations in the image data will 
naturally depend on how the image data are obtained. If the 
image data are obtained on the basis of radiographic images, 
the low frequency intensity variations may be due to, e.g., 
scattering of the electromagnetic radiation, anatomic struc- 

30 tures surrounding the illuminated part of the trabeculae such 
as cortical bone, fat tissue and muscles of varying thickness 
as well as, e.g., inhomogeneous X-ray illumination of the 
part of the bone. 
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According to the first aspect of the present invention, local 
image intensity information and variation in the local inten- 
sity are utilized to extract information relating to the 
trabecular structure of the part of the bone. In, e.g., 
5 digitized radiographic images, local image intensity informa- 
tion may be the individual pixel values, whereas variation in 
the local intensity is related to the textural information 
contained in, e.g., inhomogeneities in the image data. 

The extracted features resulting from the image manipulation 
and feature extraction procedure quantify properties of the 
trabecular structure and, thus, of the bone quality. The 
extracted features are subsequently introduced into an esti- 
mation procedure in which a predetermined relationship 
between features and bone quality enables the estimation 
procedure to estimate the bone quality of the vertebrate. 

In the present context, to "emphasize" information, such as 
magnitude information, means to give prominence to prevailing 
frequency information. This may be performed by either enhan- 
cing the prominent information or by reducing the less domi- 
20 nant information - optionally both. 

According to the first aspect of the invention, an estimate 
of the bone quality is obtained in an estimation procedure on 
the basis of a predetermined relationship between features 
obtained and reference bone quality parameters. This prede- 
25 termined relationship is typically established through stat- 
istical modelling, where explanatory variables (image fea- 
tures and optionally other explanatory features relating to 
the bone quality) are used to model corresponding reference 
bone quality data (response variable) . 

30 As described above, even though bone quality is not equal to 
bone quantity in the present context, the reference bone 
quality parameters are preferably parameters related to the 
strength of the bone rather than, e.g., BMD or BMC which, as 
explained above, relate to the density of the bone rather 
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than the strength thereof. However, to provide even higher 
diagnostic sensitivity, combination of the present invention 
with a measurement of BMC/BMD is contemplated to be useful. 
BMC/BMD measurements can be obtained either by the above - 
5 mentioned two- wavelength- technique or they may, if the image 
data are obtained on the basis of an X-ray image, be avail- 
able from the same image data, where, e.g., a standard alu- 
minum wedge has been illuminated together with the bone in 
question. 

10 At present, the most common way of obtaining non- invasive 

information about the bone quantity and bone quality is based 
on X-ray illumination. X-ray images may easily be used to 
generate the image data for use in the present invention. As 
it is preferred to perform the statistical analysis on a 

15 computer, the analogue X-ray image may be digitized and 
introduced into a computer by scanning the X-ray image. 

Naturally, the resolution of the X-ray film used to obtain 
the X-ray image will have an effect on the quality of the 
image data. Thus, it is presently preferred that the resol- 

20 ution of the X-ray film is at least 4 pairs of lines per 

centimetre, such as at least 5 pairs of lines per centimetre. 
Even though this may be sufficient for obtaining a satisfac- 
tory estimation of the bone quality, it is preferred that the 
resolution of the X-ray film is at least 10 pairs of lines 

25 per centimetre, such as at least 25 pairs of lines per centi- 
metre, preferably at least 50 pairs of lines per centimetre. 
It is possible to obtain resolutions of the X-ray film up to 
at least 100 pairs of lines per centimetre, such as at least 
250 pairs of lines per centimetre and probably as high as 500 

30 pairs of lines per centimetre, such as at least 600 pairs of 
lines per centimetre. Naturally, a resolution of this 
preferred size will ensure that a large amount of information 
is present in the X-ray image. 

Naturally, it is preferred rhat the scanning of X-ray images 
35 is performed with a sufficiently large resolution in order to 
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ensure that a minimum relevant information of the X-ray image 
is lost in the transfer to the image data. Thus, it is pre- 
ferred that the scanning has been performed at a resolution 
of at least 10 lines per cm, such as at least 25 lines per 
5 cm, preferably at least 100 lines per cm, such as at least 
200 lines per cm, such as at least 250 lines per cm. 

Furthermore, it is preferred that the resolution of the 
scanner is better than 4 true bits per pixel, such as better 
than 6 true bits per pixel, preferably better than 8 bits per 
10 pixel. 

Naturally, in order to have the scanner actually scan the 
image, it should . be able to transilluminate the radiographs. 

The ESKOSCAN 2450 from Eskofot A/S has been found to fulfil 
the above criterias and to be highly suited for use in the 
15 method of the invention. 

It is presently preferred that the background correction 
procedure comprises at least reducing or optionally removing 
low frequency information having a frequency significantly 
lower than the spacing of the projected trabeculae. As 

20 described above, this low frequent spectral content of the 
image data is typically caused by cortical bone, fat tissue 
and muscles of varying thickness as well as, e.g., 
inhomogeneous X-ray illumination of the part of the bone. 
Naturally, this type of undesirable effect is unavoidable 

25 when obtaining non- invasive image data on the basis of X-ray 
radiographs. However, other non-invasive image acquisition 
techniques, such as, e.g., MR and CT imaging may generate 
image data not automatically or not to the same degree "suf- 
fering from" this type of undesired effect. 

30 Even though it is difficult to exactly quantify the limit 
between undesired low frequency information and the desired 
higher frequency relevant information independently of the 
image resolution, it is presently preferred that information 
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having frequencies half or less than the spacing of the 
projected trabeculae is at least reduced or optionally 
removed. More preferably, information having frequencies 
being a quarter or less, such as a tenth or less than the 
5 spacing of the projected trabeculae is at least reduced or 
optionally removed. 

One preferred way of reducing or preferably removing the low 
frequency information is using a background correction pro- 
cedure comprising generating secondary image data as a result 
0 of performing a median filtering with a predetermined kernel 
size and subtracting this result from the original image 
data. One of the advantages of using a median filter is that 
this operation is edge preserving. 

Another way of reducing or preferably removing the low fre- 
5 quency information is using a background correction procedure 
comprising generating secondary image data as a result of 
performing a mean filtering with a predetermined kernel size 
and subtracting this result from the original image data. A 
mean filtering is typically much faster than the median 
0 filtering. However, the mean filtering is not edge preser- 
ving. Edges generated in this part of the process may have an 
adverse effect and generate false or erroneous information 
later in the process, leading to wrong estimates of the bone 
quality. 

5 It is typically preferred that the kernel size is at the most 
1/2 of the image data, such as at the most 1/4 of the image 
data, preferably at the most 1/10 of the image data, more 
preferably at the most 1/20 of the image data. 

A third way of reducing or preferably removing the low fre- 
0 quency information is using a background correction procedure 
comprising globally fitting a two-dimensional polynomial to 
the image data and generating background corrected image data 
on the basis of the residuals of the fitting procedure. Apart 
from potential difficulties in determining the optimal order 
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of the polynomial, which the person skilled in the art will 
know, this method potentially offers an extremely fast back- 
ground correction of image data. 

It is contemplated that a suitable order of the polynomial 
5 may be at the most 15, such as at the most 10, more preferab- 
ly at the most 5. 

Having performed a background correction of the image data, 
an image manipulation and feature extraction procedure is 
performed in order to obtain features quantifying textural 
10 properties of the trabecular structure. 

An important textural property relating to the bone quality 
is the degree of anisotropy of the projected trabecular 
pattern. This property may be described in an intuitively 
appealing way using the power spectrum of the image data. The 
15 power spectrum may be obtained using direct methods, auto- 
covariance methods or parametric methods. 

Classical spectral estimation, however, which comprises 
direct methods and auto-covariance methods, typically gives 
relatively non- cons is tent spectral estimates with a poor 
20 resolution. Alternatively, smoothing methods may be employed. 
This, however, brings about the well-known 
variance/resolution trade-off. 

Therefore, parametric spectral estimation is used according 
to the invention as this estimator is a consistent spectral 
25 estimator with superior resolution properties. 

One way of obtaining the parametric estimate of the power 
spectrum is subjecting the image data, optionally weighted 
with a window, to a Fast Fourier Transformation. This is a 
so-called direct method of which there are several variants, 
30 where the most popular probably is the so-called Welch 
method . 
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Another way of obtaining the parametric estimate of the power 
spectrum is subjecting an estimate of the auto-covariance 
function of the image data, optionally weighted with a 
window, to a Fast Fourier Transformation. 

5 As mentioned above, the parametric methods potentially offer 
a higher resolution than the classical methods. The basic 
principle of parametric methods is to identify an appropriate 
model of the image data, which is assumed to be a homogeneous 
texture, and subsequently estimate the parameters of the 
10 model. Having estimated these parameters, the parametric 
estimate of the power spectrum of the fitted model may be 
obtained. In one dimension, this type of spectral estimate 
has proven itself superior to the classical spectral 
estimates. 

15 The methods, on the basis of which the parametric estimate is 
obtained are at present preferably chosen from the group 
consisting of: causal Simultaneous Auto-Regressive Moving 
Average (SARMA) models, non- causal SARMA models, Gaussian 
Markov Random Field models and Maximum Entropy Spectral 

20 Estimates. 

When considering spatial data, such as image data, causality 
is not natural. Abandoning causality, however, leads to 
serious problems in the estimation of the parameters as the 
Least Squares estimator is no longer consistent. Therefore, 
25 it is often preferred to superimpose an artificial 

directionality on the image data. Thus, having a causal 
model, the Least Square estimation principle may be applied. 

However, a disadvantage of superimposing an artificial 
directionality on the image data is the so-called "direc- 
30 tional bias" which may be quite conspicuous in, e.g., the 
estimate of the parametric estimate of the power spectrum. 
Consequently, the use of causal models in the present context 
is not an obvious choice. This problem is partly solved using 
the -parallel -resistor" averaged spectral estimator 
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introduced by Chan in ("Two-dimensional Spectral estimation 
from Auto- regressive Models with varying Areas of Support" 
from 1981) . 

Hitherto, when considering non-causal models, one has been 
5 forced to consider the Maximum Likelihood Estimator, which 
gives a number of practical problems for the following rea- 
sons : 

The Likelihood function includes the determinant of the 
Jacobian matrix which is of the size M 2 xM 2 (where M in 
10 the 2-D situation is the side length of the image). This 

determinant is a complicated non- linear function of the 
model parameters. 

Even if assumptions are made about the image (e.g. torus) 
that offers an analytical expression for the above-men- 
15 tioned determinant, the computational load will still be 

very large. 

Close to non-stationarity, the Likelihood function is 
extremely rippled. This will often cause the non- linear 
optimization routine to find a sub-optimal solution to 
20 the estimation problem. It is not feasible/possible to 

test for this in two dimensions or more, due to the 
computational effort required. 

For the above-mentioned reasons, the Maximum Likelihood 
estimation is not a practical way of obtaining estimates of 
25 the model. Other estimation principles may, of course, be 
more practical alternatives. 

At present, it is preferred that the method on the basis of 
which the parametric spectrum is obtained, is a non- causal 
SARMA model using the so-called MORSE estimator. This estima- 
30 tor will be described further below. 

Having obtained the parametric estimate of the power spectrum 
of the image data, one or more features are extracted which 
quantify properties of the bone quality. 
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The distribution of the energy in the parametric estimate of 
the power spectrum relates directly to the anisotropy of the 
image data. Thus, several types of features may be proposed 
which quantify this anisotropy. 

5 At present, it is preferred to obtain at least one feature 
related to parameters of a contour encompassing at least a 
predetermined percentage of the energy of the parametric 
estimate of the power spectrum. 

This the contour may, e.g., is determined by 

10 (a) defining a contour, in terms of one axis in each 
dimension, around the centre of the parametric 
estimate of the power spectrum, all points on the 
contour having substantially the same distance to the 
centre of the parametric estimate of the power 

15 spectrum, and the contour encompassing less than the 

predetermined percentage of the energy of the 
parametric estimate of the power spectrum, 

(b) for each dimension of the data calculating the per- 
centage of the energy of the parametric estimate of 

20 the power spectrum encompassed by a dilated contour 

in which the axis in the dimension in question is 
increased by a predetermined distance, 

(c) increasing the axis of the contour, in the dimension 
in which the increase in energy encompassed by the 

25 dilated contour is the largest, by the predetermined 

distance and 

<d) repeating steps (b) and (c> until the percentage 
encompassed by the contour exceeds or equals the 
predetermined percentage . 



30 



The distribution of energy in the parametric estimate of the 
power spectrum is typically oblong. Tests have shown that the 
distribution of the energy in the parametric estimate of the 
power spectra of osteoporotic patients differ significantly 
from those of non-osteoporotic persons. This difference may, 
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e.g., be seen in a difference in the area covered by or the 
elongatedness of the above-mentioned contour ♦ 

One way of obtaining suitable features from the contour is to 
define the shape of the contour and subsequently derive 
5 features from parameters of the contour. 

A suitable way of defining the contour is to define the axes 
of the contour so that these are orthogonal and directed 
along principal directions of the parametric estimate of the 
power spectrum. 

10 A natural shape of this contour may be an ellipsoidal contour 
and one or more features may be derived from the semi -axes 
thereof. However, a rectangular contour may also be used. 

The above-mentioned method of defining the contour is prefer- 
ably performed so as to obtain the smallest possible contour 
15 of the desired shape. This will ensure that the contour is 
uniquely defined and that the features derived from parame- 
ters of the contour relate as well as possible to the anisot- 
ropy of the image data. 

An alternative way of extracting information in the frequency 
20 domain is to derive information from e.g. the height, width 
or total area under peaks of a smoothed periodogram. From 
this information, features relating to periodicity of the 
trabeculae may be derived. In this manner, information 
relating to a periodicity of the trabeculae is obtained by 
25 performing a frequency analysis of the image data. 

Several other ways of extracting features from the parametric 
estimate of the power spectrum may be used, such as fitting a 
Gaussian to the normalized parametric estimate of the power 
spectrum and deriving features from this Gaussian. Other 
30 methods are using rings and wedges as described by Weszka et 
al. (1976). 
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An image manipulation and feature extraction method alterna- 
tive to or supplementing the above use of the parametric 
estimate of the power spectrum is to estimate the projected 
trabecular pattern of the image data by using a Fourier 
5 method to emphasize the information in the image data 
relating to the trabecular structure and subjecting the 
manipulated image data to a morphological operation, and 
extracting features relating to the trabecular structure from 
the estimated projected trabecular pattern. 

10 The emphasizing procedure is essential in order to bring out 
the projected trabecular pattern in a robust manner. To this 
end, Fourier methods have proven suitable, and it is pre- 
ferred that the Fourier method comprises 

subjecting the image data to Fourier transformation, 
15 - performing a subsequent mathematical transformation of 
the Fourier transformed data, 

converting the transformed data back into the spatial 
domain. 

Naturally, the subsequent mathematical transformation may be 
20 linear or non-linear. 

Even though it is preferred that the subsequent mathematical 
transformation is performed by raising the magnitude of the 
Fourier transformed data to a power larger than 1, such as a 
power in the range of 1.1-10, preferably 1.5-4, such as about 
25 2.0, also other transformations, such as other types of 
filtering, may be performed. 

In order to further emphasize the information relating to the 
trabecular structure, a third mathematical transformation may 
be performed in order to preserve only part of the magnitude 
30 information of the transformed data. This may, e.g., be 

removal of at least part of the information in the magnitude 
of the spectrum. 
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Having emphasized the information in the image data relating 
to the trabecular structure, the resulting image data is 
preferably subjected to grey- scale morphological operations 
in order to prepare the data for feature extraction. 

5 At present, the purpose of the grey- scale morphological 
operations is to extract a binary representation of the 
projected trabecular pattern (PTP) . The foreground colour 
represents the trabeculae, whereas the background colour 
represents the cavities. 

10 The above PTP is preferably obtained using the so-called 
morphological top-hat operation followed by a thresholding. 
The result of these operations is preferably further 
"cleaned" by removing small isolated clusters of pixels. 

On the basis of the PTP obtained, a number of different 
15 features may be obtained. In the following, the generation of 
just a few of the vast number of possible features is 
described. 

One method of generating features from the PTP is performing 
an operation (the so-called "distance transform") comprising 
20 determining, for each background pixel, the distance accord- 
ing to a given metric from the background pixel to the 
nearest foreground pixel. The features extracted from this 
operation are indeed intuitively appealing as they relate 
directly to the inter- trabecular distance. 

25 The above metric is presently the Euclidian distance. How- 
ever, other metrics may eventually prove more suitable in the 
present context. 

Features based on the distance transform may, e.g., be 
derived from a mean value and/or standard deviation and/or 
30 the coefficient of variation and/or skewness and/or kurtosis 
of the determined distances. 
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Another method of generating features from the PTP is per- 
forming an operation (hereinafter referred to as "star area 
transform") comprising determining a measure for each back- 
ground pixel based upon determining the distance from the 
5 background pixel to the nearest foreground pixel in a number 
of given directions in the image data. 

Again, features derivable from the star area transform may, 
e.g., relate to a mean value and/or standard deviation and/or 
the coefficient of variation and/or skewness and/or kurtosis, 
10 optionally the maximum distance, of the determined measures. 

In addition to the above-mentioned features derivable from 
the image data, it may be preferred to input additional data 
into the estimation procedure. It is evident that a large 
number of features which relate to the illuminated bone, but 
15 which are not derivable from the image data, may enhance the 
precision of the estimation of the bone quality. 

As the present method may be used to estimate the bone qual- 
ity of any bone from any vertebrate, additional information 
as to the age and/or sex, and/or species, and/or race and/or 
20 the specific bone considered in the vertebrate, and/or a 
estimated Bone Mineral Density, and/or a estimated Bone 
Mineral Content, may be included in the estimation procedure. 

Even though the BMD may be introduced in the estimation 
procedure, this measure may also optionally be determined 
25 from the image data. BMD may be estimated by including data 
from a reference object in the exposure of the bone to the 
electromagnetic radiation and on the basis of the absorption 
of the electromagnetic radiation of the bone and of the 
reference object. 

30 The purpose of the estimation procedure is to estimate given 
bio -mechanical properties of the bone on the basis of the 
introduced features extracted from the image data and other 
explanatory variables. 
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One way of obtaining an estimation procedure of the above 
type is to have the estimation procedure based on a statisti- 
cal model, taking into account the correlation structure in 
the data set, so as to assign appropriate weights to the 
5 significant features in accordance with the predetermined 
relationship. 

The above-mentioned model may be determined in a number of 
ways. However, establishing a model of this type requires 
obtaining corresponding values of all relevant features and a 
10 response variable relating to the bio-mechanical property 
which is to be estimated by the estimation procedure. The 
bio-mechanical property in question may be an absolute or a 
relative measure of the bone quality. 

Examples of interesting bio-mechanical properties of bone are 
15 the mechanical bone strength, and/or a Bone Mineral Density 
measurement, and/or a Bone Mineral Content measurement, 
and/or a score value by a skilled radiologist. Naturally, the 
method of the present invention will output a value corre- 
sponding to the calibration of the method: if the method is 
20 calibrated towards a strength parameter, the output of the 
method will relate to the strength of the bone in question. 

An important purpose of the present invention is to provide a 
method for evaluating bone strength and fracture risk of 
significant clinical value, thus to provide a higher estima- 
25 tion of fracture risk than the best technique available 

today, the Double X-ray Absorptiometry Bone Mineral Content 
measurements (DXA-BMC measurements) . 

As the calibration of the estimation procedure will determine 
the estimated parameter of the bone, this parameter should be 
30 chosen in a manner so that it correlates to the bone quality. 
The assessment of the extent to which the embodiments of the 
present invention are able to estimate bone quality may be 
performed by: 
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1. Cohort -studies, in which occurrence of new fractures are 
recorded and related to the initial estimation. 



2. cross-sectional studies, in which fracture cases and 
age- and- sex -matched controls are compared to estimations 

5 obtained by the method according to the invention and option- 
ally further to DXA-BMC measurements, or 

3. by using bio-mechanical testing of bone specimens, pre- 
ceded by estimating the bone quality according to the inven- 
tion and optionally additionally DXA-BMC -measurements of the 

10 bone in question. 

Using either distal radius, lumbar spine, femoral neck, or 
any bone in which two-dimensional image data comprising 
information relating to the trabecular structure the bone can 
be obtained, coherent values for determined bone strength, 
15 DXA-BMC and the estimate according to the present invention 
may be performed on the basis of human post-mortem bone 
samples or samples from other vertebrates . 

The calibration of the method of the present invention may 
also be calibrated along the above- illustrated methods as 
20 these methods generate bone quality information and estima- 
tion on the basis of corresponding image data. 

Bone strength may, e.g., be evaluated by directly measuring 
the plasticity and maximum load of a bone in a stress- strain 
diagram (See, e.g., Lis Mosekilde) . These parameters are 
25 preferably obtained in both the direction of the vertical 
trabeculae and in the orthogonal direction in order to have 
the most complete estimate of the bone quality. 

In the method of the invention, substantially all types of 
models may be used in the estimation procedure. At present, a 
30 preferred model is selected from the group consisting of: a 
General Linear Model, a Generalized Linear Model, an Artifi- 
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cial Neural Network, a Causal Probabilistic Net or Classifi- 
cation And Regression Trees. 

As described above, the present method may be used to esti- 
mate the bone quality of bone in any vertebrate. Thus, the 
5 vertebrate may be a human, a horse, a great ape, a large ape, 
an anthropoid ape, a pig, a cow, etc, and the actual bone may 
naturally be virtually any bone in the vertebrate, such as 
radius, femur, corpus vertebrae, calcaneus, talus, os carpi, 
metatarsi, metacarpi, falanges, tibia, fibula, patella, ulna, 
10 humerus, mandible, clavicula, scapula, os coxae, os 

naviculare, os cuboideum, os cuneiform I, os cuneiform II or 
os cuneiform III. 

A second aspect of the present invention relates to a method 
for estimating the bone quality of a vertebrate, on the basis 

15 of two-dimensional image data comprising information relating 
to the trabecular structure of at least a part of a bone of 
the vertebrate, the image data being data obtained by expos- 
ing at least the part of the bone to electromagnetic radi- 
ation, the method comprising subjecting the image data to a 

20 statistical analysis comprising 



25 



a background correction procedure in which low frequent 
intensity variations not related to the trabecular struc- 
ture of the bone is reduced relative to image data 
related to the trabecular structure of the part of the 
bone, 

an image manipulation and feature extraction procedure 
comprising subjecting the image data to at least one of 
the following procedures: 

(a) obtaining an estimate of the parametric estimate 



30 



of the power spectrum of the image data and 
extracting features relating to the energy 
distribution of the parametric estimate, 



35 



(b) obtaining, on the basis of image data on which a 
Fourier method has been used to emphasize the 
information in the image data relating to the 
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trabecular structure, an estimate of a grey- level 
co-occurrence matrix and extracting at least one 
feature on the basis of the estimated co-occur- 
rence matrix, 

5 (c) obtaining an estimate of the projected trabecular 

pattern of the image data by using a Fourier 
method to emphasize the information in the image 
data relating to the trabecular structure and 
subjecting the manipulated image data to a mor- 
10 phological operation, and extracting features 

relating to the trabecular structure from the 
estimated projected trabecular pattern, 
(d) obtaining, on the basis of a frequency analysis 
of the image data, features relating to the 
15 periodicity of the trabecular structure of the 

part of the bone, 
and an estimation procedure in which the bone quality of 
the vertebrate is estimated with a Multiple Correlation 
Coefficient better than 0.5 on the basis of the derived 
20 features and optionally other features related to the 

bone or the vertebrate and a predetermined relationship 
between the features and reference bone quality parame- 
ters . 



In fact, the bone quality of the vertebrate is preferably 
25 estimated with a Multiple Correlation Coefficient better than 
0.55, such as better than 0.6. Naturally, the higher the 
Multiple Correlation Coefficient the higher the correlation 
between the actual bone quality and the estimated bone qual- 
ity. Thus, a Multiple Correlation Coefficient better than 
30 0.65, such as better than 0.7, or better than 0.8, such as 
better than 0.85 is preferred in order to have as high a 
correlation to the actual value as possible. 

In a third aspect, the present invention relates to a method 
for estimating Non-Causal Simultaneous Auto-Regressive Moving 
35 Average models in two or more dimensions, the method compris- 
ing optimizing a given direct measure of the flatness of the 
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residual spectrum of the model. In the present context 
"optimizing a given direct measure of the flatness" means 
that the flatness is maximized. 

Non- linear optimization procedures are typically used to 
maximize measures of the above type. A suitable non-linear 
procedure may comprise the following steps 

(a) generating a set of initial parameters for the model, 

(b) generating the residual spectrum of the model on the 
basis of the parameters, 

(c) obtaining the measure of the flatness of the residual 
spectrum, 

(d) obtaining a new iterate of the parameters on the 
basis of the flatness measure and a search direction 
in parameter space, 

(e) repeating steps (b) - (d) until given stop criterion is 
reached . 

A number of measures of the flatness of the residual 
parametric estimate of the power spectrum may be chosen. A 
well known measure of flatness of probability distributions 
20 is the so-called entropy. Therefore, in the present context, 
the entropy is measured on the normalized residual parametric 
estimate of the power spectrum, disregarding the DC value. 

Preferred embodiments of aspects of the invention will now be 
described with reference to the drawings wherein 

25 Pig. l shows a typical radiographic image of a wrist, 
obtained by an ESKOFOT scanner (See Example l) , 

Fig. 2 shows an extracted sub- image of the image of Fig. X, 

Fig. 3 shows the image of Fig. 2 after a median filtering 
with a 25x25 kernel size, 



10 



30 Fig. 4 shows the final background corrected image, 



WO 96/07161 PCT/DK9S/00338 

25 

Fig. 5 shows the power spectrum of the background corrected 
image of Fig. 4, 

Fig. 6 shows the restored result of the image of Fig. 4 after 
having been subjected to an emphasizing procedure, 

5 Fig. 7 shows the result of a morphological top-hat operation 
on the restored image of Fig. 6, 

Fig. 8 shows the image of Fig. 7, wherein small isolated 
groups of pixels (noise) have been removed (this image is 
referred to as the projected trabecular pattern (PTP) ) , 

10 Fig. 9 shows the result of a distance transform calculated 
for the PTP of Fig. 8, 

Fig. 10 illustrates the star area transform for a single 
pixel . 

Fig. 11 shows the result of the star area transform of all 
15 background pixels in the PTP of Fig. 8, 

Fig. 12 shows the result of the maximum distance transform of 
the PTP of Fig. 8, 

Fig. 13 shows the smallest possible ellipse containing 70% of 
the spectral energy for a 'normal' patient, 

20 Fig. 14 shows the smallest possible ellipse containing 70% of 
the spectral energy for an osteoporotic patient, 

Fig. 15 illustrates the systematic bias of the LS estimates 
for simulated isotropic textures (the vertical lines along 
the curve outline the empirical confidence intervals for the 
25 mean of the estimated parameters ( , 

Fig. 16 illustrates the systematic bias of the approximative 
ML estimates (close to non-stationarity) for simulated 
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isotropic textures (the vertical lines along the curve out- 
line the etnpirical confidence intervals for the mean of the 
estimated parameters) , 

Fig. 17 illustrates a part of the likelihood function for a 
5 first -order bilateral SAR model {given a realization of the 
process with the parameters oH3=-0 .22735 . The likelihood is 
plotted in the interval (a,/S) € [-0.275. -0.20) x t-0. 275, -0.20) ) , 

Fig. 18 illustrates that good estimates may be obtained using 
the torus -ML estimator with the approximate ML estimates as 
10 initial values combined with a stationarity check (if a 
non- stationary solution model is found, the optimizer is 
restarted with new initial values; the vertical lines along 
the curve outline the empirical confidence interval for the 
mean of the estimated parameters) , 

15 Fig. 19 illustrates that good estimates may be obtained using 
the MORSE estimator (the vertical lines along the curve is 
the empirical 95% confidence intervals for the mean of the 
estimated parameters) , 

Fig. 20 shows the objective functional of the MORSE estimator 
20 close to non- stationarity las opposed to, e.g., the ML-estim- 
ator, this objective functional is not rippled, which makes 
it much easier to perform the non- linear optimization) , 

Fig. 21 illustrates the periodogram estimator of the power 
spectrum for a normal individual (left) and for an 
25 osteoporotic individual (right) . The spectrum is shown in 

inverse so that the high intensity values are darker than the 
lower intensity values, 

Fig. 22 illustrates the parametric NSHP estimate of the power 
spectrum for a normal individual (left) and for an 
30 osteoporotic individual (right) . The spectrum is shown in 

inverse so that the high intensity values are darker than the 
lower intensity values, 
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Fig. 23 illustrates the parametric WQP estimate of the power 
spectrum for a normal individual (left) and for an 
osteoporotic individual (right) . The spectrum is shown in 
inverse so that the high intensity values are darker than the 
5 lower intensity values, 

Fig. 24 summarizes an experiment performed in order to 
describe the applicability of the MORSE estimator. 

Fig. 25 illustrates histograms based on parameters estimated 
for each of 1000 synthesized 128x128 textures, 

10 Fig. 26 illustrates optimal strengths plotted against the 
predicted strength, where age is not included in the model, 

Fig. 27 illustrates optimal strengths plotted against the 
predicted strength, where age is included in the model, and 

Fig. 2B illustrates the correlation of the logarithmic 
15 optimal fracture load and a single feature suitably chosen. 

EXAMPLE Is Textural analysis of bone structure 

In the following Example, a preferred embodiment of the 
method of the invention will be described wherein X-ray 
images of the distal radius and a vertebral body from the 

20 lumbar spine (L3) are restored using image processing. The 
restored images are used as a basis for extracting textural 
features that are shown to correlate well with the structure, 
the density and the fracture load of the trabecular bone. 
Several ways of extracting promising textural features are 

25 described. A model based on textural features (and possibly 
age and sex) is described and shown to discriminate well 
between osteoporotic and non- osteoporotic cases. 
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The X-ray images used in this Example were obtained using a 
standard X-ray instrument having a focus area of 0.6 mm 2 , a 
tube voltage of 45 kV, a focus-to-foil distance of 100 cm, a 
5 single, fine X-ray foil (it should be noted that no foil is 
required), using a double emulsion film (it should be noted 
that a single emulsion film is equally applicable) and adjus- 
ting the MAS product to obtain a suitably exposed film. 

The X-ray pictures were scanned into a computer using a 
10 ESKOFOT ESKOSCAN 2450 scanner capable of scanning up to 5000 
lines/cm. However, at present 600 lines/cm was used. This 
scanner was found highly suited for the present use. A typi- 
cal image is shown in Fig. 1. 

The demands to a scanner suitable for this use are that it 
15 should naturally be able to trans illuminate the X-ray film, 
it should preferably have a resolution better than the grains 
in the film/foil combination (so as to loose no information 
in the scanning procedure) and it should preferably generate 
grey- values of at least 8 true bits, in order to obtain a 
20 suitable dynamic range of the image data. 

Spatial Texture Models 

Preprocessing 

In this study, only the textural information is taken into 
account. Therefore a sub-image is extracted in a region of 

25 the bone containing only trabecular bone, as e.g. shown in 
Fig. 2. Low frequent intensity variations due to, e.g., 
scattering of the radiation, anatomic structures as cortical 
bone, muscles, fat tissue of varying thickness etc., are 
removed by subtracting a median filtered version of the image 

30 (shown in Fig. 3) from the image itself (a technique known as 
"unsharp filtering"). In this case a 31x31 median filter is 
used. The background- corrected image is shown in Fig. 4. This 



WO 96/07161 PCI7DK95/00338 

29 

image is very noisy and should be further processed prior to 
further transformations and subsequent feature extraction. 

The noise in the image is clearly very high- frequent, whereas 
the information which is preferably considered {the trabecul- 
5 ar structure) is of a more low-frequent nature. 

The power spectrum of the extracted, background corrected 
sub-image is shown in Fig. 5. In the following, an emphasiz- 
ing procedure will be described, which may not be generally 
applicable for removal of high-frequent noise, but tends to 
10 emphasize the dominating (in this case low-frequent) features 
in the images, while suppressing less dominating (in this 
case high frequent noise) features. This property makes it 
attractive for this specific purpose. 

Algorithm 1 The following steps outline the restoration 
15 algorithm that removes high-frequent noise and emphasize the 
relevant structures in the image. 

1. Perform a 2-D FFT and obtain the phase, X p (f r , f c ) , and 
the magnitude, X^f,., f c ) . 

2. The magnitude is raised to the power p: 

S(f r .f e ) H**<^c>l p (1> 



20 3. The histogram of S(f r , f c ) is stretched to match a Gaus- 
sian: S G (f r , f c ) 

4. A context image, C(f r , f c ) , is formed by thresholding 
S G (f r , f c ) preserving T% of the pixels. Objects less than n 
pixels are removed from the thresholded image, (note that 

25 more than T% of the energy is preserved) . 

5. A 'new' magnitude image, X N (f r , t c ) , is formed using the 
context image and the power spectrum (not the magnitude) : 
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X if / )=l 5(f " f c) i/C(/,./ e )-l (2 , 
*»ir r ,r e j | Q Q otherwise 

6. The new magnitude image X N (f r , f c ) is combined with the 
phase image X p (f r , f c ) , and the restored image, Xr(**c), is 
obtained as the real part of the inverse Fourier transform. 
It is preferred to use thresholding in the range T=l-6 and 
5 p=l-2 . 

Note that the magnitude, raised to the power p, is used to 
form the new magnitude image. This may not be obvious, but it 
has been experienced that this in general provides very 
useful results. The reason for this is, that when the 
10 magnitude is raised to a given power larger than one, the 

parts of the spectrum that contain most of the energy tend to 
be emphasized whereas spectral components with less energy 
tend to be weakened further. The result of the emphasizing 
procedure using T=2 and p=2 is shown in Fig. 6. 

15 

Extraction of the PTP 

In the following, various image transformations and textural 
features, extracted from transformed images, will be dis- 
cussed. Several features are based upon an image subsequently 
20 referred to as the projected trabecular pattern (PTP) which 
is obtained as outlined below: 

Algorithm 2 The steps outlined below describe the desired 
transformations of the restored image in order to measure a 
relevant set of features on the resulting PTP. 

25 1. Perform the grey-scale morphological top-hat opera- 

tion on the restored image. The top-hat operation is 
the original image minus an opened version of the 
image. This image is subsequently thresholded. The 
result of these operations is shewn in Fig. 7. 
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2. This image clearly contains a lot of noise-pixels. 
Small isolated objects (in the present Example all 
objects less than 10 pixels) are removed. The result- 
ing PTP is shown in Fig. 8. 

5 Extraction of Statistical Features 

Naturally, it is preferred to extract features that, in some 
sense, quantify the properties {such as structure and 
density) of the trabecular structure. Not surprisingly, the 
meshes of the PTP-net are much larger, and often broken, for 
10 osteopenic individuals than for non-osteoporotic individuals. 
This difference may be quantified in several ways. Such 
properties are, e.g., reflected in the average distance from 
each background pixel to the nearest foreground pixel 
(bone-pixel) in the projected trabecular pattern. 

15 The size of the dark meshes in the pattern is a relevant 

parameter that is conveniently quantified by calculating the 
so-called distance transform for the projected trabecular 
pattern. The result of the distance transform on the PTP is 
shown in Fig. 9 and is denoted xo( r ' c > • 

20 In addition, for each background pixel in the projected 
trabecular pattern, a line is drawn from the pixel, in a 
given direction, until a foreground pixel { 'bone- pixel ' ) is 
found. The distance between these pixels is measured. This is 
done for a fixed number of directions as illustrated in Fig. 

25 10. After this, the average distance which is referred to as 
the star area, denoted x A (r,c), or the maximum star length 
which is denoted xx( r » c) be calculated. These images are 
shown in Fig. 11 and Fig. 12 respectively. 

A number of statistical moments may be calculated from the 
30 resulting images from the measured average, standard 

deviation, coefficient of variation and skewness for all 
background pixels: 
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A. = 



(3) 



d 



(4) 




(5) 



V <X.(r,c)-M 5 



(6) 



where the subscript denotes different transforms considered: 

• D: Distance transform 

• M: Mean value of star length 

• S: Standard deviation of star length 
5 • X: Maximum star length 

• .A: Star area 

We briefly outline the interpretation of the considered 
parameters : 

The (l. and <f. features are either related to the general mesh 
10 size (for D, M, A and X) or the degradation of the mesh (S) . 
For example, larger mesh sizes are reflected in a larger 
value of {l . . Thus, for increasing values of JT. and , the 
texture index should decrease. 

The cv. feature is somewhat different. It measures the width 
15 of the empirical distribution in the transformed images. If 
all the distances are large, cv. would be small. The same 
thing would be the case, if all the distances were small. 
Typically, however, we find that the empirical distribution 
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is wide, only if the PTP mesh is broken. Thus, for increasing 
values of cv., the texture index should typically decrease. 

The 7. feature is also related to the mesh size. If the 
histogram of distances is heavier in the right side, y.<0. If 
the empirical distribution is symmetrical, y*0 . For large 
meshes in the PTP, the distribution of distances would 
typically be skewed for large values. Thus, for increasing 
values of 7., the texture index should decrease. 

Tvorturai Futures E xtr act ed in the sp atUi frequency Domain 

One would intuitively think that features measured in the 
frequency domain should carry relevant information about the 
trabecular structure. The trabecular structure for strongly 
osteoporotic patients is clearly more low frequent than that 
for 'normal' patients. In the present section, the considered 
Fourier features are not measured on the projected trabecular 
pattern. Instead, the features are measured on the spectrum 
of the background corrected image. 

Preprocessing 

Before estimating the spectrum, non-stationaries of the 
extracted image are removed by subtracting a 31x31 median 
filtered version of the image from itself as described above. 

Spectral Estimation 

As the main interest is directed toward the relatively low- 
frequent properties of the power spectrum, a weighed estimate 
of the dispersion may be used, where pMu.v) is the 
normalized power spectrum taking a weight function into 
account : 
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p'lu.v).— p(u ' v) - ^'"'^ ... 



and 
w(u, v) i 



(8) 



Given a positive definite symmetrical matrix 




the eigenvalues are: 
. A°l+°Wl*l-*l) 2 +**tt {x0) 

1 o 



2 



(11) 



5 and the corresponding (normalized) eigenvectors 



Pi- 



(12) 



An alternative estimate of the power spectrum may be obtained 
from the following process: 

The probably best known estimator of the power spectrum is 
the periodogram estimator: 
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M-lM-l 
Af n>-0 n-0 ' 




(13) 



where (u,v) are spatial frequencies. 

It is well known that this estimate is a non- consistent 
estimate of the power spectrum. For example, when the samples 
from a Gaussian white noise process are used, the periodogram 
5 estimate has a variance that is proportional to the variance 
of the white noise process. 

One approach to variance reduction is to segment the data, 
such as in non -overlapping blocks, where the spectrum is 
estimated in each block, using this periodogram estimator, 
10 and obtaining the final estimate by averaging the 

periodograms of the blocks. In this manner, the variance is 
decreased by a factor equal to the number of blocks. Of 
course, this is at the expense of frequency resolution. 

In the following, an iterative process is suggested for 
15 estimating the dispersion matrix, which estimate is more 

sensitive to the anisotropics reflected in the low frequency 
parts of the orientation of the texture. 

Algorithm 3 The dispersion matrix of the approximating 
Gaussian is estimated iteratively in the following way: 



20 



1. E 0 = I 



2. i - 1 

3 . repeat 

3a. Estimate the dispersion: 




(14) 



where 
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(15) 



u v 



fit»J2£ v 2 p(u, v) v(u, V) 



(16) 



U V 



d uv = ]£ 5") UV P ( u, v) *r( u, v> 




(17) 



and 



w(u. v) - (u, v) 




(18) 



3b. i = i-1 
until | | E A - Ei.Jl.wx < 6 

From both these techniques, the principal direction of the 
5 power spectrum may be found using: 



This principal direction may be used to find the ellipse with 
the smallest possible area, containing a fixed fraction of 
the energy of the original (non- weighted) power spectrum. The 
ellipse, rotated with the angle or around the origin, is 
10 expressed in polar coordinates: 




(19) 



RW 2 - 



a 2 b 2 



(20) 



b 2 cos 2 (6-a) +a 2 sin 2 (9-a) 



where a and b are the semi -major and semi -minor axis, 
respectively, and where the above-mentioned \'s are the 
eigenvectors . 



In the present Example, the semi -major and semi -minor axes 
15 may be found using the following scheme: 
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Algorithm 4 



The steps outlined below are used to estimate 



the semi -major and semi-minor axis in the ellipse with the 
smallest possible area, containing a given fraction of the 
spectral energy. The angle with which this ellipse should be 
5 rotated around the origin is found as outlined above. 

1. Start with the semi -minor and semi -major axis suffi- 
ciently small, e.g., a = b = 10. 



2. Increase a and b with one, and calculate the fraction 
of energy in each ellipse (rotated with the angle or 



3. The ellipse with the largest fraction of energy is 
chosen . 

4. Repeat the steps 2 and 3 until the requested fraction 
of energy is reached. 

The features considered here are the estimated values: a, b, 
a-b, b/a as well as the features \ x , X 2 , X 2 /X 1( P=a uv /<r u <x v , 
o v 2 /a u 2 , and: 



10 



around the origin) . 




(21) 



(22) 



v /(o 2 u -ot) 2 -4ot 



uv 



<U 2 *V a )p<!J, V) 



(23) 



u v 



-££p(u, v) log <p(u,v)) 



(24) 



u v 



In Fig. 13, the estimated ellipse containing 70% of the 
20 energy {and the direction) is shown for a 'normal' patient. 
In Pig. 14, the 70% ellipse is shown for an osteoporotic 
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patient. Clearly, the frequency content for the osteoporotic 
patient is a lot more low- frequent . This means that the 
energy is concentrated near the centre of the image and, 
consequently, the areas of the ellipses are very different 
5 although they contain the same fraction of energy. 

Sparta,} Jnceraccion tfQteXs 

Above, it is described that a powerful set of features may be 
extracted in the spatial frequency domain. However, the 
performance of features extracted in the frequency domain 
10 depend heavily on the spectral estimator used. The 

periodogram is a non- consistent spectral estimator. This 
problem may be overcome by segmenting the data in non- 
overlapping segments, estimating the periodogram in each 
segment and finally averaging over these periodograms . 

15 The variance reduction is directly proportional to the number 
of segments. However, this is at the expense of frequency 
resolution. 

In the following, we will consider two types of causal 
Simultaneous AutoRegressive (SAR) models to obtain parametric 
20 spectral estimates; the Non- Symmetrical Halfplane (NSHP) 
model and the Quarter- Plane (QP) model. 

Assume a wide sense stationary, zero mean, stochastic process 
defined on a rectangular grid of pixel sites Q = {s«(s r ,s c ): 
0*s r sM, 0ss c sN} , obeying a SARMA model defined as: 



25 where r=(r r ,r c )€ or N^. and are the regions of 
support for the parameter arrays {<t>} and (0), respectively. 
We also refer to these sets a neighborsets which, in some 
sense, defines the order of the model. 




(25) 
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The innovation process £ s is assumed to be an IID N(0,Z f 2 ) or 
Gaussian white noise sequence. If = 0, we have a pure SAR 
model, and if ~ 0, we have a pure SMA model. 

In the following, we will focus on pure autoregressions , for 
5 the following reasons: The use of mixed SARMA models would 
probably give more parsimonious models, but the models 
required to capture the nature of the considered texture are 
of a such low order that parsimony in not really an important 
issue. It also turns out that the important feature presently 
10 is the shape of the spectrum in terms of e.g. peak locations.* 
Also, the LS estimator for the pure SAR model is a set of 
linear equations, whereas non- linear optimization is required 
in order to obtain estimates of mixed SARMA models. 

The (non- normalized) spectral density function is given as 
15 the squared modulus of the frequency response function 
multiplied by the noise variance: 



,2 rc^g 



Py(f)=0l 



* I 1 * E ©x*xp<-2iciT T -f)f 



(26) 



f2 (i+flTr(*W f )N^(;W/) 2 



where the spatial frequencies (f r ,f c ) = fcQ f , and 



(27) 



The vector notation is defined as: 



(28) 



(29) 
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C(N. ) f =£Ql[cos (to (f , r) ) , re*/. ] (30) 
S(N. ) t =ool [sin (Q (f,r) ) ,r6N.] (31) 
q ( f . r) =2* tr r -f e +r c '£ c ) (32) 



So far, the question of causality has not been addressed. In 
two dimensions, the notion of causality is not natural. 
Abandoning causality, however, leads to a number of 
difficulties, one of them being that the Least Squares 
5 estimator of the parameters is no longer consistent. Hence, 
it has become popular to impose an artificial directionality 
on the image, thus obtaining a causal model (in fact, 
recursively computatable is a more correct term) . 



The two types of causal models typically considered in the 
10 literature are the Quarter Plane support models (QP) and the 
Non- Symmetrical Half Plane support models (NSHP) . In fig. 21, 
parametric estimates using the NSHP model are illustrated. It 
is seen that the spectrum is slightly elongated in a certain 
direction, when compared to the corresponding periodogram 
15 estimate seen in Fig. 22. This phenomenon is due to 
directional bias which is an effect that is even more 
pronounced for the QP mode. In order to solve this problem, a 
Weighted QP (WQP) estimator has been proposed: 



P AV (f r ,f e )> 



(33) 



where the subscripts and -+ simply refer to the considered 
20 quadrant, i.e. the first and the second quadrant. This 

parametric estimator is widely considered to be among the 
best parametrical spectral estimators, but it suffers from 
the drawback that it has no model interpretation. 

In fig. 23, a WPQ estimate is illustrated in which no 
25 directional bias is seen. 
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Using these estimates, features of the type described above, 
that is: \ 1( X 2 , XjAi. p=ff rc /a r ff c , <r c 2 /o T 2 , and: 



— = ) 



(34) 



(35) 



;(o 2 f -a 2 c ) 2 -4a^ 
may be used. 

Also, as cr describes the principal orientation of the 
5 approximating Gaussian, this direction may be used to find 
the above-mentioned energy ellipses. Thus, also the semi- 
major and semi-minor axes of these ellipses a and 6, 
respectively as well as the features a-b, b/a may be used. 

F Kr ranrAnn of features relying %o r?erx<?<tiGitY 

10 An alternative method of deriving features from the 

preprocessed image of Fig. 6 is to again subject this image 
to the two dimensional Fourier transformation and extract the 
Fourier power spectrum (the periodogram spectrum) . 

It is presently preferred that the resolution of the image is 
15 reduced using a Gaussian pyramid, thus, preserving as much of 
the information as possible in the high resolution prior to 
the Fourier transform. 

Feature Extraction 



It is not surprising that the power spectrum comprises peaks 
20 corresponding to periodicity of the vertical trabeculae. It 
is, however, fortunate that also peaks - however small - 
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corresponding to a periodicity of .the horizontal trabeculae 
may be discerned from the power spectrum. 

Information relating to the density of trabeculae, their 
thickness and the distance between them may be extracted from 
5 the height of these peaks relating to the periodicities, the 
total area under these peaks and/or the "width" thereof. 

At present, it is preferred not to extract the peak positions 
from the Fourier power spectrum, .as the periodogram or the 
power spectrum are very noisy due to the fact that the 
10 periodogram estimator of the power spectrum is non- 

consistent. Instead, a combined algorithm is preferred, where 
the parametric spectrum is used to identify these peak 
positions, and the smoothed periodogram is used to estimate 
the associated energy etc. 

15 In the power spectrum, a total of four peaks are typically 
present - or more precisely, only two peaks, as the 
parametric spectrum is symmetric. A line is drawn through the 
two largest peaks (which relate to the vertical trabeculae) 
through the DC-value. From the middle of. this line, a line is 

20 drawn perpendicularly thereto also through the DC- value. It 
is inherent from the model that the peaks of the horizontal 
trabeculae are positioned perpendicularly to those of the 
vertical trabeculae. The maxima on the new line will 
correspond to the peaks of the horizontal trabeculae. 

25 This method is preferred in order to be independent and 
insensitive to any rotation of the image. 

Smoothing of the spectrum may be performed using a kxk 
kernel, such as a 5x5 kernel using, e.g., a mean filter, a 
median filter, a Gaussian filter etc. 

30 From this spectrum, the peak heights, the volume under the 
peaks a.s.o. may be determined and used in the prediction of 
the bone quality. 
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In addition, it is contemplated that also the curvature of 
the peaks in the parametric power spectrum may be used as a 
predictor of bone quality, as this curvature also relates to 
deviations of the periodicity. 

5 EXAMPLE 2: Study of feature performance 

A data set was used to investigate different image features 
and compare their performance. Radiographic images were 
obtained from 97 individuals and each image was ranked by 
experienced radiologists. The analysis of this data set is 
10 described in the following. 

Each patient was assigned the score +1 if the radiologist 
considered the patient to be 'normal' and the score -1 in 
case of accelerated bone loss (osteoporotic) . The four radio- 
logists agreed upon the diagnosis in 52 cases, and only these 

15 are used in this preliminary study. In this study of feature 
performance, a general linear model with the score as 
response variable and age, sex and textural features as 
explanatory variables is used. Furthermore, rigorous statis- 
tical tests of the significance of the individual features is 

20 performed using a backwards elimination procedure. 

The features used have been described above, and are all 
based on first, second and third order statistics measured on 
transforms of the projected trabecular pattern (PTP) or 
fitted ellipses on the Fourier power spectrum of the restored 
25 images. 

The fitted models are of the type 

Si - P 0 ♦ M^i * M e *i * Mi.i * * € i < 38> 

where s i is the score of the i'th patient, age i the age of 
the i'th patient, eex i the sex of the i'th patient (+1: Male, 
-1: Female in order to make both sexes contribute equally to 
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the mean) and fj #i the j ' th image feature of the i'th 
patient. It is assumed that € is i.i.d. gaussian white noise. 
Under this assumption, the LS estimates are equal to the ML 
estimates. Test statistics are calculated for each parameter 
5 in the model. Using a backwards elimination approach, the 
least significant parameter of the model is eliminated, and 
the model is subsequently re -estimated with the remaining 
parameters. This is continued until no more parameters can be 
eliminated from the model. 

10 Features ba sed on the projected trabecular pattern 

Out of 12 image features, 9 were eliminated using the back- 
wards elimination procedure. The order in which they were 
eliminated is summarized in the table below: 



Feature 
number 


Feature 
name 


1 


■2 




4 b 




7 




9 


1* i /i£) 1 


















2 


°D 








| • 






























4 


to 






• 












h 


Us 










• 








6* 


<7.C 




















cv s 


• 
















8 


7s 


















9 














• 






10 


vsi 
















• 


11 


cv M 




• 














12 


7M ! 








1 









15 The boldfaced features marked with an asterisk are signifi- 
cant parameters that could not be eliminated. 

The resulting model is given by 

Si " Po * + M e *i * Pl^.i + P4^0.i * Ms + € J (39) 



A summary of the resulting model is given below: 
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Call: lm( formula 



score (idx4) - 

age + sex[idx4] + f[idx4 f 1) + f [idx4, 3] 
+ f[idx4, 6)) 



Residuals : 



5 Min 


1Q 


Median 


3Q 


Max 


-3.755 


-0.8549 


-0.06046 


0.7199 


3.854 


Coefficients 












Value 


Std. Error 


t value 


Pr(>|t|) 


(Intercept) 


50.2650 


12.2399 


4.1067 


0.0002 


10 age 


-0.0995 


0.0141 


-7.0427 


0.0000 


sex[idx4) 


1.1937 


0.2579 


4.6281 


0.0000 


f[idx4, l] 


-11.7739 


2.3699 


-4.9681 


0.0000 


f[idx4, 3] 


-37.2419 


11.0728 


-3.3634 


0.0016 


f[idx4, 6] 


2.7775 


1.2370 


2.2454 


0.0297 



15 Residual standard error: 1.746 on 45 degrees of freedom 
Multiple R-Squared: 0.7516 



20 



25 



Correlation of Coefficients: 
(Intercept) age 
age 0.0063 
sex[idx4j -0.1329 
f(idx4, 1) -0.6991 
f[idx4, 3] -0.9837 
f[idx4, 6] -0.0287 



sexlidx4] f|idx4.11 £11(1x4.3) 



0.2430 
-0.1573 -0.0231 
-0.0538 0.1220 0.7305 
0.1205 0.1072 -0.6269 -0.0918 



Most of the explanatory effect is contained in the age and 
the sex. The rest of the variability {approximately 16% 
point) is explained by the textural features. It may be 
desired to transform some of the features, or even derive new 
ones . 



From the above, it may be concluded that the considered 
30 statistics measured on at least the distance transform and 
the star area of the projected trabecular pattern, are quite 
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strong features. This investigation, however, does not sug- 
gest the maximum distance transform as an optimal feature. 

Four4gr Fee tvreg 

In the following, four features extracted from the power 
5 spectrum will be considered. The only feature eliminated from 
this set of features is the length of the semi-major axis 
(a) . This is in good agreement with what would be expected, 
as the erosion of the trabecular structure in osteoporotic 
patients is mainly in the vertical direction. The performance 
10 of the resulting model is comparable to the results obtained 
using features based on the distance transform. 

Again the data were analyzed in the statistical software 
package Spins, and the following output was obtained for the 
resulting model: 

15 Residuals: 

Min 1Q Median 3Q Max 





-2.852 - 


1.114 -0. 


.3153 0.6817 


4 


.854 






Coefficients : 
















Value 


Std. Error 


t value 


Pr(>|t|) 


20 


(Intercept) 


8.1977 


2.9152 


2 


.8120 


0.0073 




ages [idx4 , 2] 


-0.1086 


0.0153 


-7 


.1036 


0.0000 




sex ( idx4 ] 


1.2674 


0.2711 


4 


.6747 


0.0000 




freq[idx4, 2] 


8.6433 


2.S788 


3 


.3516 


0.0016 




freq[idx4,3] 


-0.2102 


0.0660 


-3 


.1845 


0.0026 


25 


freq{idx4,4] 


-68.6948 


24.0649 


-3 


.6857 


0.0006 



Residual standard error: 1.843 on 45 degrees of freedom 
Multiple R-Squared: 0.7232 
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Correlation of Coefficients: 









>u(kU4) 




ages [idx4,2] 


0.0165 








sex(idx4) 


0.0103 


0.2150 






freq[idx4,2] 


-0.2471 


-0.1272 


0.1435 




freq(idx4,3) 


0.2169 


0.1303 


-0.1491 


-0.9972 


freq{idx4,4] 


0.1145 


0.086b 


-0.1543 


-0 .9880 



0.9845 



Thus, from the results of the above it may be seen that a 
Multiple Correlation Coefficient better than 0.7 may be 

10 obtained using only either the features extracted from the 
power spectrum or the features extracted from the PTP; only 
part of the features offered by the method of the invention. 
Thus, it is contemplated that even higher Multiple Correla- 
tion Coefficient may be expected when these groups of fea- 

15 tures are combined in the estimation procedure. In this 

connection it should be mentioned that Multiple Correlation 
Coefficients of this size are quite satisfactory. 



EXAMPLE 3: The MORSE Estimator for Mixed Non-causal SARMA 
Models 

20 The concept of texture is widely recognized as being of great 
importance in many fields of application, such as, e.g., 
medical imaging, industrial quality inspection and remote 
sensing. Visual properties of a given texture are, however, 
not very tangible and many attempts has been made to extract 

25 good and meaningful features in order to describe such prop- 
erties. 

One such approach is to model the spatial correlation struc- 
ture using texture models. The two mainstream types of models 
are the Gaussian Markov Random Field (GMRF) models and the 
30 Simultaneous Auto- Regressive Moving Average (SARMA) models. 
The latter type of models is an attempt to extend the very 
powerful ARMA models (known from time-series analysis and 1-D 
signal processing) to two dimensions. 
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Whereas l-D signals are inherently causal, this is not a 
natural constraint in two dimensions. Abandoning causality, 
however, introduces severe problems in important modelling 
steps as identification and estimation. Therefore, it is very 
5 popular to impose an artificial directionality on data and 
thereby obtaining consistent estimates by minimizing the sum 
of squared estimation errors - i.e. Least Squares estimates. 

The analysis of non- causal SARMA models has been severely 
constrained, mainly due to the fact that the non-consistency 

10 of the LS estimator forces one to consider a very complicated 
Likelihood function. As a consequence, the analysis has been 
concentrated solely on pure SAR models. The toroidal ML 
estimator and an approximative ML estimator has been formu- 
lated for pure SAR models. It can be shown that the perform- 

15 ance of the approximative estimator may not be adequate when 
the considered textures is close to non-stationarity . 

Considering the exact ML estimates for toroidal SAR models 
will, in general, give very good results. It is shown, how- 
ever, that the likelihood function is extremely rippled in 

20 the vicinity of non-stationarity. This will, in many cases, 
cause the optimization routine to get 'caught' in a sub- 
optimal solution. This problem may indeed be serious as 
visual textures are often characterized by a strong positive 
correlation structure yielding estimates close to non-statio- 

25 narity. 

Another problem, which is related to the Maximum Likelihood 
estimator and the traditional Maximum Entropy estimators, is 
that these tend to overemphasize the high-frequency content 
of the considered. 

30 Thus, a new, entropy based, estimator is introduced. Instead 
of maximizing the entropy of the associated model spectrum 
(as in MEM spectral estimation introduced in "A new algorithm 
for two-dimensional maximum entropy power spectrum estima- 
tion" Lim, J. S. & Malik, N. A. (1981) IEEE Transactions on 
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Acoustics, Speech and Signal Processing, 29, 401-413), the 
entropy of the residual spectrum is maximized. This estimator 
is named the MORSE estimator (Maximum Of Residual Spectral 
Entropy) . This estimator is well behaved, even close to 
5 non-stationarity, which makes it much simpler to perform the 
non- linear optimization of the involved objective functional 
(w.r.t. the parameters). It is possible to write the objec- 
tive functional for mixed SARMA models (i.e. including Moving 
Average parts) and obtain very good estimates. To the best of 
10 our knowledge there are no results in the literature regard- 
ing inference in non -causal SARMA models. 

The maximization of the residual spectral entropy seems, 
intuitively, to be a reasonable criterion as the underlying 
assumption of the simultaneous models is that the innovation 

15 process is white. It is obvious that the optimality of a 
given model is closely related to the chosen criterion of 
optimality, and the choice of entropy is by no means the only 
one to be imagined. Note that no assumptions are made 
concerning the driving noise (innovation process ????) 

20 sequence. We only require it to have a symmetrical 
distribution. 

yhP non -causal SARMA model 

Consider a wide sense stationary, zero mean bilateral SARMA 
process sampled on a square MxM grid of pixels s»(r,c)efl 



ys+E^-r " C S * £ 6 ' € *-r (40) 

25 where r»(r 1( r r ) £ N. and N. defines the order of the model. 
The only assumption that is made about the innovation process 
(y.) is that it is independent identically distributed random 
variables with a symmetrical density. 



WO 96/07161 



50 



PCT/DK95/00JJ8 



Contrary to the GMRF models, the bilateral SARMA model does 
not have to be symmetrical, but opposite neighbours must be 
associated with identical parameters {i.e. tf> r » #. r . This is 
necessary to ensure the identif iability of the parameters. 
5 For the sake of simplicity, the following is restricted to 
consider only symmetrical models . 

In order to facilitate the notation when it is desired to be 
specific about the order of the model, the notation outlined 
in the following definitions is introduced: 

10 Definition 1 A rectangle with side lengths 2a+l and 2b+l, 
rotated by an angle ot around the origin is denoted R(a,b,a) . 
Now the following is defined 

• R(a,b) = R(a, b,0) 

• R(a) * Ria, a, 0) 

Definition 2 An ellipse with semi -major axis a and semi- 
minor axis b, rotated by an angle a around the origin is 
15 denoted E(a ( b,or) . The following is defined 

• E(a,i>)= E(a,b,0) 

• E(a) ■ E(a,a,0) 

Thus, the order of a given model may be defined in terms of 
these geometrical shapes as, e.g., SARMA (E (l) , E (l) ) which 
means that the auto-regressive part, as well as the moving 
average part, is the four nearest neighbours. 



20 



If the bilateral model is considered, the stationarity condi- 
tion is greatly weakened; all that is required is that no 
poles of the transfer function falls on the unit circle: 
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(l + EWrV'oj (43) 



for all |z r | =1 and |z c | = 1. 

The transfer function is given as the 2-D z- transform of the 
impulse response function. 

H{2 - ) *c> _ ^IlU)^ 6 ^ 2 ^ (44 ) 



and the frequency response function is obtained as the 2-D 
5 Fourier transform of the impulse response function 



where 

(46) 



is the spatial frequency and 



£<N.) f( „ = £Ql(cos<2irf (s) T r) ,r€tf.] 



In time- series analysis and signal processing, the successful 
use of ARMA models for spectral estimation has promoted an 
10 interest in the extension of these methods to two dimensions. 
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As for ARMA models, the 2-D parametric spectrum is propor- 
tional to the squared magnitude of the frequency response: 



S(f(s))=ol\'H(f(s))\ 2 (48) 

where o € 2 is the variance of the innovation sequence. 

Using this approach, the problem of spectrum estimation 
5 becomes one of fitting an appropriate model to the observed 
data. Mostly, parametric spectral estimation in 2-D has been 
attempted using unilateral models. This typically leads to 
directional bias. By using a bilateral model, this particular 
problem should be solved, provided a practically applicable 
10 estimation principle is provided. Our preliminary results 
show that the MORSE estimation principle is a very useful 
tool for obtaining parametric spectral estimates without 
directional bias . 

Different approaches have been made in order to obtain esti- 
15 mates for the non-causal SAR model. There have, to the best 
of our knowledge, been no results regarding estimation of 
mixed non-causal SARMA models. Even for the pure SAR model 
the estimation of the parameters is not a simple problem. 
Below, different estimators for the non- causal SAR model will 
20 be illustrated and their shortcomings discussed. 

The Non- causal (or bilateral) Simultaneous Auto- Regress ions 
(SAR) models, sometimes referred to as NCSAR, were originally 
introduced by P. Whittle in "On stationary processes in the 
plane". Biometrika, 41, 343-449. The problem of non-consist- 

25 ency of the LS estimator is addressed in this paper, and a 
large sample approach based on spectral methods was pres- 
ented. First, this method will be of use only for very simple 
models and secondly, many practical situations involve small 
samples for which the discounting of edge effects may not be 

30 negligible. In order to illustrate how the non- consistency of 
the LS estimator affects the estimates, a simple experiment 
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is carried out. For an isotropic nearest neighbour model, 100 
realizations are generated for different parameters. In each 
case, the parameters are estimated (as if the model was 
anisotropic). In Fig. 15, the average of the estimated para- 
5 meter is plotted versus the 'true' parameter {i.e. the para- 
meter with which the texture was simulated) . It is seen that 
the LS estimator does not provide satisfactory estimates. In 
fact, the order of magnitude is off by a factor two in many 
cases . 

10 Kashyap, R. & Chellappa, R. proposed an approximation to the 
exact (toroidal) ML- estimator in "Estimation and choice of 
neighbours in spatial interaction models of image", IEEE 
Transactions on Information Theory, 29(1), 60-72. The approx- 
imation is based on a Taylor- series expansion of the deter- 

15 minant of the Jacobian (which becomes block circulant with 

circulant blocks under the torus assumption) . Only first- and 
second- order parts (in the parameters) are included. The 
improvement of the results is definitely vast, as shown in 
Fig. 16 where the above experiment is repeated. Nevertheless, 

20 as the bias is very pronounced close to non-stationarity, it 
is seen that this estimator may not be adequate in many 
practical situations where the textures tend to lie close to 
non-stationarity, and it is highly preferred to obtain a more 
suitable estimator not having these disadvantages. 

25 Grunkin investigated the exact (toroidal) ML-estimator in "On 
The Analysis of Image Data Using Simultaneous Interaction 
Models" in IMSOR, ph.d thesis # 67, Institute of Mathematical 
Statistics and Operations Research, Technical University of 
Denmark, Lyngby, 223pp. The likelihood function involved is 

30 complicated and computationally very expensive. Also, the 
likelihood function is extremely rippled in the vicinity of 
non-stationarity as shown in Fig. 17. Whether the non-linear 
optimization routine is 'caught' in a sub-optimal solution or 
not is very much dependent upon the initial guess and the 

35 parameters of the optimization algorithm (such as, e.g., step 
length). Thus, it may take several re -estimations to be 
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convinced that the correct solution is found. This, of 
course, is merely a technical problem, and in the simple 
experiment described above, these problems can be overcome 
and the exact (toroidal) ML- estimates yield very good results 
5 as shown in Pig. 18. 

As mentioned above, the non- consistency of the LS estimator 
is leads to non-sensical results. The practical application 
of the exact toroidal ML estimator will, in many cases, be 

10 limited, as the likelihood function tends to be extremely 
rippled in the vicinity of non-stationarity . Thus, it may 
require several re -estimations to ensure that the 'true' 
optimum is found. Also, as pointed out earlier, the ML esti- 
mator may overemphasize the high- frequent content of the 

15 texture. In addition, the computational load involved in 
obtaining estimates for mixed SARMA models using the ML 
estimation principle is prohibitive for all practical 
purposes . 

The approach taken here is to consider the spectral content 
20 of the residuals. Finally, it will probably be extremely 
difficult to obtain ML estimates of mixed non- causal SARMA 
models . 

In this Example, a novel estimator is described that gives 
very good estimation results, has a response surface that {in 
25 our experience) is very well behaved and allows for inference 
in mixed SARMA models. 

As the underlying assumption is that the residuals should be 
white, i.e. uncorrelated, all frequencies should be represen- 
ted, with the same weight, in the spectrum of the residual 
30 process. If the residual spectrum is normalized, the entropy 
of the residual spectrum obtains its maximum when the spec- 
trum is as flat as possible. This observation is the basis of 
the Maximum Of Residual Spectral Entropy (MORSE) estimator. 
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Consider again the difference equation defining the system: 



y.*£**y« * ^ + E 0 rS-x (49 > 



The output can be viewed as the convolution of the innovation 
process with the impulse response of the system. Thus, in the 
frequency domain, the following representation is found: 

y(f(s))= tf(f(s))- Elf is)) <50) 



5 where Y(f(s)) is the Fourier transform of the output, E(f{s)l 
is the Fourier transform of the (non- observable) innovation 
process and M(f (s) ) is the Fourier transform of the impulse 
response function, which is also known of the frequency 
response function. 

10 Using equation (50) , the power spectrum of the innovation 
process may be expressed in terms of the power spectrum of 
the observations and the parametric spectrum (assuming that 
the parameters are known). Thus it is found that: 



\Elf(s) ) P. I y(f (S) » I' -|r<f <«) ) |» (jggg^ «d (51) 



A normalization function depending upon the observations and 
15 the parameters may be defined: 



The parameters are subsequently found by maximizing the 
entropy of the residual spectrum. 
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Definition 3 The Maximum Of Residual Spectral Entropy 
(MORSE) estimates are found by maximizing the entropy func- 
tion M¥,£aR'£ma> with respect to the parameters fi^ and fi^: 



=arg max £.(^,§^,9^) {53) 



where 

(54) 



An estimate of the residual variance is obtained as 
dl =N(x. fi^, fi^) /M 2 (55) 



Clearly, this function is strongly non- linear in the parame- 
ters. The optimization may be carried out using, e.g., a 
Quasi -Newton method. In this work, the IMSL implementation of 
10 the BFGS method is applied. 

SxperitnentaX results 

As illustrated above, the exact (torus-) ML- estimates give 
very good results. The problem here is of a more technical 
nature. Firstly, the involved likelihood function (or objec- 
15 tive functional) is extremely complicated and computationally 
expensive to evaluate. Secondly, the likelihood function is 
very rippled in the vicinity of non-stationarity . 
Consequently any optimization- routine might easily 'get 
caught' in a sub -optimal solution. 

20 Above, important problems as centrality, consistency, 

efficiency of the proposed estimator have not been addressed. 
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Obviously, questions as these are difficult to answer 
analytically due to the untraceable form of the involved 
expressions. Below, these questions will be addressed 
empirically through simulation and estimation experiments. 

5 Estimation of pure SAR models 

In Fig. 19 , the objective functional of the MORSE estimator 
is shown in the same region as the ML-estimator . It is clear 
that the objective functional of the MORSE estimator is much 
smoother - i.e. there are no ripples. It is our experience 
10 that the objective functional of the MORSE estimator never 
exhibits ripples in parameter space - not even outside the 
admissible region. This makes the non-linear optimization a 
rather straight -forward task that can be undertaken 
automatically by most standard software packages. 

15 The experiment carried out for the various estimators above 
in connection with Fig. 17 is also carried out for the MORSE 
estimator. The result is shown in Fig. 20, and it is seen 
that the results are very acceptable. Note that the variance 
decreases as non-stationarity. is approached. This is due to 

20 the fact that the poles in the transfer function are 

approached. Such poles (or near-poles) create sharp peaks in 
the objective functional and thereby a much smaller variance. 

Estimation of mixed SARMA models 

A major advantage of this approach is that an estimator for 
25 mixed SARMA models is obtained 'for free'. In order to illus- 
trate that good results are obtained for the mixed SARMA 
models as well, the following experiment is conducted. A set 
of parameters are chosen for a SARMA (E( 1.0) , E (1, o) ) model. 
Furthermore, these parameters are chosen not too far away 
30 from non-stationarity so that strong positive spatial auto- 
covariances are obtained. Now 100 simulations are carried out 
and the histogram of each parameter is plotted with the 
•true' parameter marked as shown in Fig. 20. 
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Another test is one where 1000 realizations of size 64x64 
pixels and 1000 realizations of size 128x128 pixels are 
generated - also not too far from non-stationarity . For each 
of these realizations, the MORSE estimator was used to obtain 
5 an estimate of the model parameters. The empirical 95%- 
confidence intervals are calculated for the mean of each 
parameter. In fig. 24, it is seen that the known values of 
the parameters fall nicely within these confidence intervals 
further indicating that the estimator is central. It should 
10 also be noted that the width of the confidence intervals 
decrease as the size of the synthesized texture samples 
increase. This indicates that the estimator is consistent. 

Thus, it is seen that very satisfactory results are obtained 
also for the mixed non-causal SARMA model. 

15 Distributional properties of the estimator 

It might intuitively be expected that the distribution of the 
estimated parameters be Gaussian if the driving noise was 
Gaussian. As this is hard to prove formally, this is made 
probable through experiments. Histograms of the parameters, 

20 estimated from the 1000 128x128 simulated textures described 
above, are shown in Fig. 25. The hypothesis that the shown 
histograms represent samples of Gaussian estimates will now 
be formally tested using the x 2 test. For each parameter, the 
estimates are divided into K classes, where K is selected as 

25 K = 1+3 .31og 10 (N) (where N is the number of observations). In 
order to make the test as strong as possible, the classes are 
selected such that the expected number of observations is 
equal in each class. This is done by generating the splits as 
the (i-100/k)%-quantile (i=l K-l) . The test statistic is 

30 now computed as: 
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Where ObSi is the number of observations in the i'th class. 
This test statistic is distributed as x 2 {K-3) (since both 
mean and standard deviation should be estimated) . The tests 
are summarized in the below table, where it is seen that the 
5 distributional hypothesis is easily accepted for all 
parameters . 





&i 


mean of 0; 


sdev of 9; 


Z 


P(Z > X 2 (8)} 


AR(-l t 0) 


-0.22314 


-0.22308 


0.00757 


5.44 


0.709 


A*(0,-1) 


•0.12118 


-0.12126 


0.00808 


6.94 


0.543 


MA(-l,Q) 


0.11535 


0.11521 


0.00772 


7.512 


0.483 


AfA(O.-l) 


0.24066 


0.24083 


0.00759 


6.038 


0.643 



The principle on which this estimator lies is intuitively 
10 very appealing. The empirical results obtained for this 

estimator is, indeed, very promising. Using this estimator, 
also non- causal SARMA models may be considered. 

EXAMPLE 4: Test of the presently preferred method 

In the present example, a test was performed in order to 
15 evaluate the applicability of the present method in the 

prediction of bone strength on the basis of X-ray pictures of 
the bone. 

Pete ffgt 

The data set consisted of 20 samples of the third lumbar 
20 vertebrae harvested from cadaveric specimens each of which 

were exposed to X-ray irradiation in order to obtain an X-ray 
picture and a BMD measurement of each bone. 

Subsequently, an anterior and a posterior core sample of 
trabecular bone was extracted from each bone, and the optimal 
25 fracture load was determined for each of these samples using 
a three-point test setup as known gsx s£. 
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The X-ray pictures were scanned with a spatial resolution of 
600 lines per cm. In the final resolution, an area is 
selected centrally positioned between the positions of the 
drilled-out samples. The positions of the core samples may be 
5 partly or fully comprised in the selected area. 

It was decided to choose the mean optimal fracture load of 
the two core samples of a bone as the strength of the bone. 
This largely corresponds to the strength of the bone in the 
central part as determined by linear interpolation of the 
10 strengths of the two samples. The results of the present test 
seems to confirm that this decision is reasonable. 

Pre-nrocessina of the imaae material 

In order to have as good images as possible in the final, 
much smaller resolution, the scanned X-ray images had a large 
15 spatial resolution (600 lines per centimetre) . 

The present images were so large that a scale reduction using 
a Burt-Adelson Gaussian pyramid was impossible due to the 
computer having too little memory. Thus,, as a preliminary 
reduction was performed by firstly smoothing the image using 
20 a 3x3 mean filter and subsequently discarding every other 
pixel in the image. 

After this preliminary reduction, the image had a size 
reducible by the Burt-Adelson Gaussian pyramid to obtain the 
final resolution (150 lines per centimetre) . 

25 Subsequent to the reduction of the image, the image was 

background corrected by subtracting a 31x31 median filtered 
version of itself. 

From the corrected image, a Region Of Interest (ROD having a 
size of 256x256 pixels was selected. Preferably, this ROI was 
30 selected close to the center of the bone due to the fact that 
a large gradient of the thickness of the trabeculae was 
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1. Search through the list of possible features until the 
first significant feature is identified. 

2. Continue the search until another significant feature is 
identified and included in the model. 

5 3. Test whether previous features should be eliminated. 

Steps 2 . and 3 . are repeated until all features have been 
investigated. Naturally, a level of significance should be 
decided for the described tests. In the present 
investigation, this level is chosen at p=0.10, as the amount 
10 of data does not make a very conservative test possible. 

It became clear that the optimal fracture load related 
exponentially to the parameters of the model, whereby the 
model was performed as a regression with the logarithm to the 
optimal fracture load as a dependent variable. 

15 Choosing the parameters according to the above principle, 
only three significant parameters were found: 3 (See 
Algorithm 4), Formula (23) squared, and n s 2 (See Formula 
(3)) . 

In the following, these parameters are evaluated in two 
20 models, one using age as a variable and one not using age. 

Mnrifillina without acre 

A General Linear Model not using age as a descriptive 
variable was made. This model will, thus, seek to predict the 
optimal fracture load on the basis of the extracted 
25 structural indices. The result is summarized in the below 
program printout: 

Residuals : 

Min 1Q Median 3Q Max 

-0.4327 -0.2409 0.03446 0.1625 0.5578 
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Coefficients : 

Value 

11.2930 
-70.8473 

5 Ms' -0.0287 
(Formula 23) 2 143625.7515 



(Intercept) 
a 

2 



Std. Error t value Pr(>|t|) 

1.8815 6.0020 0.0000 

21.1912 -3.3432 0.0044 

0.0037 -7.8621 0.0000 

53274.0853 2.6960 0.0166 



Residual standard error: 0.2901 on 15 degrees of freedom 
Multiple R-Squared: 0.8158 

F-statistic: 22.14 on 3 and 15 degrees of freedom, the p- 
10 value is 9.17e-06 



Correlation of Coefficients: 

(Intercept) a n s 2 
a -0.9873 

Ms 2 -0.5412 0.4485 

15 (Formula 23) 2 0.4347 -0.5455 -0.1318 



Explaining model : 



F^«exp( -70. 8473* (Formtfia23) 2 -0. 0287-5 (57) 
♦143625. 7 515^1*11. 2930) 



No. 


log(Msrd) 


log (Predicted) 


Rciidnii 


238 


4.574778 


4.302492 


0 .272286112 


240 


4.259244 


4.168912 


0.090332350 


242 


4.516842 


4.742086 


.0.225244419 


247 


4.677491 


4.800312 


•0.122620941 


250 


3.508197 


3.473737 


0.034459077 


254 


4.070737 


4.402031 


•0.331293705 


258 


2.501764 


2.864725 


-0.362961146 


259 


3.724428 


3.682320 


0.042107301 


265 


4.443404 


3.885623 


0.557761393 


267 


5.382806 


5.098572 


0.284233332 


269 


4.099266 


3.949321 


0.149945096 


279 


4.434097 


4.125558 


0.308539021 


280 


3.844675 


4.277384 


.0.432709063 


281 


3.372283 


3.310585 


0.061698025 


289 


3.823039 


3.806975 


0.014063294 


293 


3.779999 


3.604959 


0.175040294 


296 


3.527713 


3.764313 


.0.256600239 


298 


3.896777 


3.886981 


0.009796051 


313 


4.397531 


4.666183 


.0.268651834 
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In Fig. 26 a plot of the logarithmic optimal fracture load 
against the predicted strength is illustrated. In this 
figure, the measurement points are represented by their 
number in the database. 



5 At this point it is, naturally, interesting to know whether 
the BMD measurements performed on all samples may take part 
as a statistically significant feature contributing to 
increasing the multiple R- squared. The result may be seen 
from the below printout 



10 Residuals: 

Min 1Q Median 3Q Max 

-0.4342 -0.2389 0.0385 0.1615 0.5598 



Coefficients : 
15 (Intercept) 

a 

Formula 23 
squared 
20 BMD 



Value Std. Error t value Pr(>|t|) 

11.3316 2.2940 4.9397 0.0002 

-71.1532 23.9423 -2.9719 0.0101 

-0.0288 0.0044 -6.5255 0.0000 

144161.1908 57644.4790 2.5009 0.0254 

-0.0190 0.5953 -0.0319 0.9750 



Residual standard error: 0.3002 on 14 degrees of freedom 
Multiple R- Squared: 0.8158 

F- statistic: 15.5 on 4 and 14 degrees of freedom, the p-value 
is 8.43e-05 



25 Correlation of Coefficients: 

(Intercept) a /i s 2 Formula 23 

squared 

a -0.9798 
li s 2 -0.6661 0.5587 

30 Formula 23 0.5070 -0.5949 -0.2583 

squared 

BMD -0.5285 0.4009 0.5154 -0.2914 
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It is obvious that BMD is not a significant feature and that 
it can be eliminated from the model. Using BMD as the only 
descriptive parameter, a multiple R-squared of R 2 «0.2142 is 
obtained indicating that BMD is a poor indicator of the 
5 quality of trabecular bone. Thus, the aspects of fracture 
risk relating to trabecular bone are poorly described by BMD. 



WfelXinq with flge 



Using age as a describing parameter, the following model is 
obtained: 

10 Residuals: 

Min 1Q Median 3Q Max 

-0.3777 -0.1385 0.02544 0.1517 0.328 



15 



20 



Coefficients: 
(Intercept) 

a 

Ms 2 

Formula 23 

squared 

AGE 



Value Std. Error t value Pr(>|t|) 

11.8882 1.6616 7.1545 0.0000 

-66.3085 18.5326 -3.6858 0.0024 

-0.0236 0.0039 -6.1152 0.0000 

152951.1278 46677.8586 3.2767 0.0055 

-0.0160 0.0067 -2.3827 0.0319 



Residual standard error: 0.2533 on 14 degrees of freedom 
Multiple R-Squared: 0.8689 

P-statistic: 23.2 on 4 and 14 degrees of freedom, the p-value 
is 4.707e-06 

H S 2 Formula 23 

squared 
I 

-0.0617 

-0.5613 -0.0838 



25 Correlation of Coefficients: 
(Intercept) a 



a -0.9658 

ji s 2 -0.3585 

30 Formula 23 0.4408 
squared 

AGE -0.1503 



0.4028 
-0.5378 

-0.0575 
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i^«exp ( -68 . 3085* (Formula21 ) 2 -0 . 0236 -<3 

-152951. 127 8ii|-0. 016 -AGE+11 . 8882) 



No. || log (Hard) 


log (Predicted) 


Reiidual 


238 


4.574778 


4.293879 


0.280898703 


240 


4.259244 


4.262127 


-0.002882869 


242 


4.516842 


4.505306 


0.011535494 


247 


4.677491 


4.766867 


.0.089375665 


250 


3.508197 


3.413527 


0.094669316 


254 


4.070737 


4.435918 


-0.365180333 


258 


2.501764 


2.852028 


.0.350263751 


259 


3.724428 


3.625144 


0.099283320 


265 


4.443404 


4.238432 


0.204971765 


267 


5.382806 


5.273592 


0.109213882 


269 


4.099266 


3.964857 


0.134408820 


279 


4.434097 


4.106100 


0J27997388 


280 


3.844675 


4.032239 


•0.187563938 


281 


3.372283 


3.394572 


.0.022288673 


289 


3.823039 


3.797599 


0.025439958 


293 


3.779999 


3.610914 


0.169084738 


296 


3.527713 


3.829923 


-0.302210802 


298 


3.896777 


3.656788 


0.239969248 


313 


4.397531 


4.775257 


.0.377726598 



In Fig. 27 a plot of the logarithmic optimal fracture load 
against the predicted strength is illustrated. In this 
5 figure, the measurement points are represented by their 
number in the database. 

Discussion 

Prom the result of the present example, it is clear that the 
optimal fracture load of trabecular bone may be predicted 

10 using only textural parameters - and even with a high 

precision. Three image features are used: Formula 23 squared, 
/i s 2 , and & describing the density of the horizontal 
trabeculae (the thickest) , the break down of the trabecular 
structure, and the global trabecular density, respectively. 

15 In addition, it is seen that neither sex nor BMD can add 
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information to the description obtained using the above three 
image features. 

However, it is seen that age is a significant parameter. 
Including age in the model increases the multiple R- squared 
5 to R 2 =0.8689. It should be noted that the two extreme 
observations (258 and 267) contribute quite a lot to the 
value of R 2 . Omitting these observations gives R 2 «0.7538, 
including age. which should still be regarded as a high 
multiple R- squared. 

10 However, as is described in connection with Algorithm 1, the 
parameters used in restoring the image will a large effect on 
the applicability of the features extracted from the 
resulting image. In fact, using p=1.4 and T*4 instead of 
p=T=2 in algorithm 1, the feature c H alone gives a R 2 of 

15 0.92. 



In Fig. 28, the correlation between <x M and the logarithm to 
the optimal fracture load is illustrated for a test based on 
the picture and optimal fracture load - material described 
above. Even if the extreme point in the lower right corner of 
20 the figure is removed, a correlation of 0.88 is obtained. 

Thus, a very large difference may be seen depending on the 
actual manner of obtaining the PTP of the image. 
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CLAIMS 



1. A method for estimating the bone quality of a vertebrate, 
on the basis of two-dimensional image data comprising infor- 
mation relating to the trabecular structure of at least a 
5 part of a bone of the vertebrate, the image data being data 
obtained by exposing at least the part of the bone to elec- 
tromagnetic radiation, the method comprising subjecting the 
image data to a statistical analysis comprising 



a background correction procedure in which low frequent 
10 intensity variations not related to the trabecular 

structure of the bone is reduced relative to image data 
related to the trabecular structure of the part of the 
bone, 

an image manipulation and feature extraction procedure 
15 wherein at least the local image intensity information 

as well as variation in the local intensity are utilized 
to extract information related to the trabecular 
structure of the part of the bone, the image 
manipulation and feature extraction procedure comprising 
20 subjecting the image data to at least one of the 

following procedures: 

(a) obtaining an estimate of the parametric estimate 
of the power spectrum of the image data and 
extracting features relating to the energy 

25 distribution of the parametric estimate, 

(b) obtaining, on the basis of image data on which a 
Fourier method has been used to emphasize the 
information in the image data relating to the 
trabecular structure, an estimate of a grey- level 

3 0 co-occurrence matrix and extracting at least one 

feature on the basis of the estimated co- 
occurrence matrix, 

(c) obtaining an estimate of the projected trabecular 
pattern of the image data by using a Fourier 

35 method to emphasize the information in the image 

data relating to the trabecular structure and 
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subjecting the manipulated image data to a 
morphological operation, and extracting features 
relating to the trabecular structure from the 
estimated projected trabecular pattern, 
5 (d) obtaining, on the basis of a frequency analysis 

of the image data, features relating to the 
periodicity of the trabecular structure of the 
part of the bone, 
and an estimation procedure in which the bone quality of 
10 the vertebrate is estimated on the basis of the derived 

features and optionally other features related to the 
bone or the vertebrate and a predetermined relationship 
between the features and reference bone quality parame- 
ters . 

15 2. A method according to claim 1, wherein the electromagnetic 
radiation is X-ray radiation. 

3. A method according to any of the preceding claims, wherein 
the image data is scanned from an X-ray film. 

4. A method according to claim 3, wherein the resolution of 
20 the X-ray film is at least 4 pairs of lines per centimetre, 

such as at least 5 pairs of lines per centimetre. 

5. A method according to claim 4, wherein the resolution of 
the X-ray film is at least 10 pairs of lines per centimetre, 
such as at least 25 pairs of lines per centimetre, preferably 

25 at least 50 pairs of lines per centimetre. 

6. A method according to claim 5, wherein the resolution of 
the X-ray film is at least 100 pairs of lines per centimetre, 
such as at least 250 pairs of lines per centimetre, preferab- 
ly at least 500 pairs of lines per centimetre, such as at 

30 least 600 pairs of lines per centimetre. 

7. A method according to claim 3, wherein the scanning has 
been performed at a resolution of at least 10 lines per cm. 
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8. A method according to claim 7 , wherein the resolution is 
at least 25 lines per cm. 

9. A method according to claim 8, wherein the resolution is 
at least 100 lines per cm, such as at least 200 lines per cm, 

5 such as at least 250 lines per cm. 

10. A method according to any of the preceding claims, where- 
in the background correction procedure comprises at least 
reducing or optionally removing low frequency information 
having a frequency significantly lower than the spacing of 

10 the projected trabeculae. 

11. A method according to claim 10, wherein information 
having frequencies half or less than the spacing of the 
projected trabeculae is at least reduced or optionally 
removed . 

15 12 . A method according to claim 11, wherein information 

having frequencies being a quarter or less than the spacing 
of the projected trabeculae is at least reduced or optionally 
removed . 

13. A method according to claim 12, wherein information 

20 having frequencies being a tenth or less the spacing of the 
projected trabeculae is at least reduced or optionally 
removed . 

14. A method according to any of claims 10-13, wherein the 
background correction procedure comprises generating second- 

25 ary image data as a result of performing a median filtering 
with a predetermined kernel size and subtracting this result 
from the original image data. 

15. A method according to any of claims 10-13, wherein the 
background correction procedure comprises generating second - 

30 ary image data as a result of performing a mean filtering 
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with a predetermined kernel size and subtracting this result 
from the original image data. 

16. A method according to any of claims 10-13, wherein the 
background correction procedure comprises globally fitting a 

5 two-dimensional polynomial to the image data and generating 
background corrected image data on the basis of the residuals 
of the fitting procedure. 

17. A method according to claim 16, wherein the order of the 
polynomial is at the most 15, such as at the most 10, prefer- ' 

10 ably at the most 5. 

18. A method according to claim 14 or 15, wherein the kernel 
size is at the most 1/2 of the image data, such as at the 
most 1/4 of the image data, preferably at the most 1/10 of 
the image data, more preferably at the most 1/20 of the image 

15 data. 

19. A method according to any of the preceding claims, where- 
in the parametric estimate of the power spectrum is estimated 
by using a method selected from the group consisting of 
direct methods, auto-covariance methods and parametric 

20 methods. 

20. A method according to claim 19, wherein the image data, 
optionally weighted with a window, is subjected to a Fast 
Fourier Transformation in order to generate the parametric 
estimate. 

25 21. A method according to claim 19, wherein the auto- 
covariance function estimated on the basis of the image data, 
optionally weighted with a window, is subjected to a Fast 
Fourier Transformation in order to generate the parametric 
estimate . 



30 22. A method according to claim 19, wherein one of the 
methods on the basis of which the parametric estimate is 
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obtained is chosen from the group consisting of: causal 
Simultaneous Auto -Regressive Moving Average (SARMA) models, 
non-causal SARMA models, Gaussian Markov Random Field models 
and Maximum Entropy Spectral Estimates and wherein the order 
5 of the model is identified, the parameters of the model are 
estimated and the spectrum of the estimated model is gener- 
ated. 

23. A method according to claim 22, wherein the method on the 
basis of which the parametric estimate is obtained, is a non- 
10 causal SARMA model using an estimator according to claim 50. 

24. A method according to any of the preceding claims, where- 
in at least one feature is related to parameters of a contour 
encompassing at least a predetermined percentage of the 
energy of the power spectrum. 

25. A method according to claim 24, wherein the contour is 
determined by 

(a) defining a contour, in terms of one axes in each 
dimension, around the centre of the power spectrum, all 
points on the contour having substantially the same 
distance to the centre of the power spectrum, and the 
contour encompassing less than the predetermined 
percentage of the energy of the power spectrum, 

(b) for each dimension of the data calculating the per- 
centage of the energy of the power spectrum encompassed 
by a dilated contour in which the axis in the dimension 
in question is increased by a predetermined distance, 

(c) increasing the axis of the contour, in the dimension in 
which the increase in energy encompassed by the dilated 
contour is the largest, by the predetermined distance 
and 

(d) repeating steps (b) and (c) until the percentage 
encompassed by the contour exceeds or equals the 
predetermined percentage . 



20 



25 
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26. A method according to claim 24 or 25, wherein the contour 
is ellipsoidal, optionally spherical. 

27. A method according to claim 25 or 26, wherein the axis of 
the contour are orthogonal and are defined along principal 

5 directions of the power spectrum. 

28. A method according to claim 1, wherein procedure (c) is 
used and wherein the Fourier method comprises 

subjecting the image data to Fourier transformation, 
performing a subsequent mathematical transformation of 
10 the Fourier transformed data, 

converting the transformed data back, into the spatial 
domain . 

29. A method according to claim 28, wherein the subsequent 
mathematical transformation is linear or non- linear. 

15 30. A method according to claim 29, wherein the subsequent 
mathematical transformation is performed by raising the data 
to a power larger than 1, such as a power in the range of 
1.1-10, preferably 1.5-4, such as about 2.0. 

31. A method according to any of claims 28-30, wherein a 
20 third mathematical transformation is performed in order to 

preserve only part of the magnitude information of the trans- 
formed data. 

32. A method according to any of claims 28-31, wherein the 
emphasized information in the image data relating to the 

25 trabecular structure is subjected to a grey- scale morphologi- 
cal operation. 

33. A method according to claim 32, wherein the grey-scale 
morphological operation is a top-hat operation wherein all 
pixels fitting in the top-hat are assigned a value different 

30 from a specified background pixel value. 
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34. A method according to claim 33, wherein the image data 
which has been subjected to the top-hat operation is thres- 
holded so as to assign a foreground pixel value to the pixels 
having a pixel value different from the background pixel 

5 value. 

35. A method according to any of claims 32-34, wherein the 
image data which has been subjected to the grey- scale mor- 
phological operation is subjected to further noise reduction 
to substantially remove isolated single pixels and pixel 

10 groups smaller than a predetermined threshold number of 
pixels having a pixel value different from the background 
pixel value. 

36. A method according to any of claims 32-35, wherein an 
operation is performed comprising determining, for each 

15 background pixel, the distance according to a given metric 
from the background pixel to the nearest foreground pixel, 
and wherein at least one feature is extracted from the result 
of the operation. 

37. A method according to claim 36, wherein at least one 
20 feature relates to a mean value and/or standard deviation 

and/or the coefficient of variation and/or skewness and/or 
kurtosis of the determined distances from all background 
pixels to the nearest foreground pixel. 

38. A method according to any of claims 32-35, wherein an 
25 operation is performed comprising determining a measure for 

each background pixel based upon determining the distance 
from the background pixel in a number of given directions in 
the image data to the nearest foreground pixel and at least 
one feature is extracted from the result of the operation. 

30 39. A method according to claim 38, wherein at least one 
feature relates to a mean value and/or standard deviation 
and/or the coefficient of variation and/or skewness and/or 
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kurtosis of the determined measures for all background 
pixels . 

40. A method according to claim 38 or 39, wherein the measure 
for each background pixel is a mean value and/or standard 

5 deviation and/or the coefficient of variation and/or skewness 
and/or kurtosis and/or maximum distance of the determined 
distances in the given directions. 

41. A method according to any of the preceding claims, 
wherein at least one feature related to the bone or to the 

10 vertebrate, such as age and/or sex, and/or species, and/or 
race and/or the specific bone considered in the vertebrate, 
and/or an estimated Bone Mineral Density of the bone, and/or 
an estimated Bone Mineral Content of the bone, is included in 
the estimation procedure. 

IS 42. A method according to claim 41, wherein the Bone Mineral 
Density is estimated by including data from a reference 
object in the exposure of the bone to the electromagnetic 
radiation and where the Bone Mineral Density is estimated by 
correlating the absorption of the electromagnetic radiation 

20 of the bone and of the reference object. 

43. A method according to any of the preceding claims, 
wherein the estimation procedure is based on a statistical 
model, taking into account the correlation structure in the 
data set, so as to assign appropriate weights to the 
25 significant features in accordance with the predetermined 
relationship . 

44 A method according to any of the preceding claims, 
wherein the reference bone quality is an absolute or relative 
bone quality of bones for which image data have been 
30 processed in accordance with any of the preceding claims. 

45. A method according to claim 44, wherein the reference 
bone quality is determined on the basis of measurement of the 
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mechanical bone strength, and/or Bone Mineral Density 
measurement, and/or on the basis of a Bone Mineral Content 
measurement, and/or on the basis of a score of the bone 
strength by a skilled radiologist. 

5 46. A method according to any of the preceding claims, 

wherein the predetermined relationship between the features 
and reference bone quality parameters is established on the 
basis of reference bone quality parameters and features 
extracted according to any of the preceding claims. 

10 47. A method according to claim 46, wherein the predetermined 
relationship is defined in terms of a model taken from the 
group consisting of: a General Linear Model, a Generalized 
Linear Model, an Artificial Neural Network, a Causal 
Probabilistic Net and Classification And Regression Trees. 

15 48. A method according to any of the preceding claims, 
wherein the vertebrate is a human. 



49. A method according to any of the preceding claims, 
wherein the vertebrate is a horse. 



50. A method according to any of the preceding claims, 
20 wherein the vertebrate is a large ape, a great ape or an 
anthropoid ape. 



51. A method according to 
wherein the vertebrate is 

52 . A method according to 
25 wherein the vertebrate is 



any of the preceding claims, 
a pig. 

any of the preceding claims, 
a cow. 



53. A method according to any of the preceding claims, 
wherein the bone is taken from the group consisting of 
radius, femur, corpus vertebrae (LI, L2, L3, L4, L5, Tl, T2, 
T3, T4, T5, T6, T7, T8, T9, T10, Til, T12, CI, C2, C3, C4, 
30 C5, C6, C7) , calcaneus, talus, os carpi, metatars, metacarpi, 



WO 96/07161 PCT/DK95/00338 

77 

falanges, tibia, fibula, patella, ulna, humerus, mandible, 
clavicula, scapula, os coxae, os naviculare, os cuboideum, os 
cuneiform I, os cuneiform II, os cuneiform III, os sacrum, os 
coccygis. 

5 54. A method for estimating the bone quality of a vertebrate, 
on the basis of two-dimensional image data comprising 
information relating to the trabecular structure of at least 
a part of a bone of the vertebrate, the image data being data 
obtained by exposing at least the part of the bone to 
10 electromagnetic radiation, the method comprising subjecting 
the image data to a statistical analysis comprising 

a background correction procedure in which low frequent 
intensity variations not related to the trabecular 
structure of the bone is reduced relative to image data 
15 related to the trabecular structure of the part of the 

bone, 

an image manipulation and feature extraction procedure 
comprising subjecting the image data to at least one of 
the following procedures: 
20 (a) obtaining an estimate of the parametric estimate 

of the power spectrum of the image data and 
extracting features relating to the energy 
distribution of the parametric estimate, 

(b) obtaining, on the basis of image data on which a 
25 Fourier method has been used to emphasize the 

information in the image data relating to the 
trabecular structure, an estimate of a grey- level 
co-occurrence matrix and extracting at least one 
feature on the basis of the estimated co- 
30 occurrence matrix, 

(c) obtaining an estimate of the projected trabecular 
pattern of the image data by using a Fourier 
method to emphasize the information in the image 
data relating to the trabecular structure and 

35 subjecting the manipulated image data to a 

morphological operation, and extracting features 



WO 96/07161 



78 



PCT/DK9S/00338 



relating to the trabecular structure from the 

estimated projected trabecular pattern, 
(d) obtaining, on the basis of a frequency analysis 

of the image data, features relating to the 
5 periodicity of the trabecular structure of the 

part of the bone, 
and an estimation procedure in which the bone quality of 
the vertebrate is estimated with a Multiple Correlation 
Coefficient better than 0.5 on the basis of the derived 
10 features and optionally other features related to the 

bone or the vertebrate and a predetermined relationship 
between the features and reference bone quality 
parameters . 

55. A method according to claim 54, wherein the bone quality 
15 of the vertebrate is estimated with a Multiple Correlation 

Coefficient better than 0.55, such as better than 0.6, 
preferably better than 0.65, such as better than 0.7, 
preferably better than 0.8, such as better than 0.85. 

56. A method for estimating Non-Causal Simultaneous Auto- 
20 Regressive Moving Average models in two or more dimensions, 

the method comprising optimizing a given direct measure of 
the flatness of the residual spectrum of the model . 



57. A method according to claim 56, wherein the optimization 
of the given flatness measure is obtained by 
25 (a) generating a set of initial parameters for the model, 

(b) generating the residual spectrum of the model on the 
basis of the parameters, 

(c) obtaining the measure of the flatness of the residual 
spectrum, 

30 (d) obtaining a new iterate of the parameters on the basis 
of the flatness measure and a search direction in 
parameter space, 
(e) repeating steps (b)-{d) until given stop criterion is 
reached . 
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58. A method according to claims 56 or 57, wherein the 
measure of flatness of the residual spectrum is obtained as 
the entropy of the normalized residual power spectrum, dis- 
regarding the DC value. 
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AMENDED CUIUS 

[received by the International Bureau on 13 February 1996 (13.02.96); 
original claims 1-58 replaced by 
amended claims 1-60 02 pages)] 
1. A method for estimating the bone quality of a vertebrate, 
on the basis of two-dimensional image data comprising infor- 
mation relating to the trabecular structure of at least a 
part of a bone of the vertebrate, the image data being data 
obtained by exposing at least the part of the bone to elec- 
tromagnetic radiation, the method comprising subjecting the 
image data to a statistical analysis comprising 



a background correction procedure in which low frequent 
intensity variations not related to the trabecular struc- 
ture of the bone is reduced relative to image data 
related to the trabecular structure of the part of the 
bone, 

an image manipulation and feature extraction procedure 
wherein at least the local image intensity information as 
well as variation in the local intensity are utilized to 
extract information related to the trabecular structure 
of the part of the bone, the image manipulation and fea- 
ture extraction procedure comprising subjecting the image 
data to at least one of the following procedures: 

(a) obtaining an estimate of the parametric estimate 
of the power spectrum of the image data and 
extracting features relating to the energy 
distribution of the parametric estimate, 

(b) obtaining, on the basis of image data on which a 
Fourier method has been used to emphasize the 
information in the image data relating to the 
trabecular structure, an estimate of a grey- level 
co-occurrence matrix and extracting at least one 
feature on the basis of the estimated co- 
occurrence matrix, 

(c) obtaining an estimate of the projected trabecular 
pattern of the image data by using a Fourier 
method to emphasize the information in the image 
data relating to the trabecular structure and 
subjecting the manipulated image data to a 
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morphological operation, and extracting features 
relating to the trabecular structure from the 
estimated projected trabecular pattern, 
(d) obtaining, on the basis of a frequency analysis 
5 of the power spectrum of the image data, a 

smoothed periodogram of the image data, or an 
estimate of the parametric estimate of the power 
spectrum of the image data, features relating to 
the periodicity of the trabecular structure of 
10 the part of the bone, 

and an estimation procedure in which the bone quality of 
the vertebrate is estimated on the basis of the derived 
features and optionally other features related to the 
bone or the vertebrate and a predetermined relationship 
15 between the features and reference bone quality parame- 

ters. 

2. A method according to claim 1, wherein the electromagnetic 
radiation is X-ray radiation. 

3. A method according to any of the preceding claims, wherein 
20 the image data is scanned from an X-ray film. 

4. A method according to claim 3, wherein the resolution of 
the X-ray film is at least 4 pairs of lines per centimetre, 
such as at least 5 pairs of lines per centimetre. 

5. A method according to claim 4, wherein the resolution of 
25 the X-ray film is at least 10 pairs of lines per centimetre, 

such as at least 25 pairs of lines per centimetre, preferably 
at least 50 pairs of lines per centimetre. 

6. A method according to claim 5, wherein the resolution of 
the X-ray film is at least 100 pairs of lines per centimetre, 

30 such as at least 250 pairs of lines per centimetre, preferab- 
ly at least 500 pairs of lines per centimetre, such as at 
least 600 pairs of lines psr centimetre. 
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7. A method according to claim 3, wherein the scanning has 
been performed at a resolution of at least 10 lines per cm. 

8. A method according to claim 7, wherein the resolution is 
at least 25 lines per cm. 

5 9. A method according to claim 8, wherein the resolution is 
at least 100 lines per cm, such as at least 200 lines per cm, 
such as at least 250 lines per cm. 

10. A method according to any of the preceding claims, where- 
in the background correction procedure comprises at least 

10 reducing or optionally removing low frequency information 
having a frequency significantly lower than the spacing of 
the projected trabeculae. 

11. A method according to claim 10, wherein information 
having frequencies half or less than the spacing of the 

15 projected trabeculae is at least reduced or optionally 
removed . 

12. A method according to claim 11, wherein information 
having frequencies being a quarter or less than the spacing 
of the projected trabeculae is at least reduced or optionally 

20 removed. 

13. A method according to claim 12, wherein information 
having frequencies being a tenth or less the spacing of the 
projected trabeculae is at least reduced or optionally 
removed. 

25 14. A method according to any of claims 10-13, wherein the 
background correction procedure comprises generating second- 
ary image c* .'.a as a result of performing a median filtering 
with a pr -ermined kernel size and subtracting this result 
from the ginal image data. 
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15. A method according to any of claims 10-13, wherein the 
background correction procedure comprises generating second- 
ary image data as a result of performing a mean filtering 
with a predetermined kernel size and subtracting this result 

5 from the original image data. 

16. A method according to any of claims 10-13, wherein the 
background correction procedure comprises globally fitting a 
two-dimensional polynomial to the image data and generating 
background corrected image data on the basis of the residuals 

10 of the fitting procedure. 

17. A method according to claim 16, wherein the order of the 
polynomial is at the most 15, such as at the most 10, prefer- 
ably at the most 5. 

18. A method according to claim 14 or 15, wherein the kernel 
15 size is at the most 1/2 of the image data, such as at the 

most 1/4 of the image data, preferably at the most 1/10 of 
the image data, more preferably at the most 1/20 of the image 
data. 

19. A method according to any of the preceding claims, where - 
20 in the parametric estimate of the power spectrum is estimated 

by using a method selected from the group consisting of 
direct methods, auto-covariance methods and parametric 
methods. 

20. A method according to claim 19, wherein the image data, 
25 optionally weighted with a window, is subjected to a Fast 

Fourier Transformation in order to generate the parametric 
estimate. 

21. A method according to claim 19, wherein the auto- 
covariance function estimated on the basis of the image data, 

30 optionally weighted with a window, is subjected to a Fast 
Fourier Transformation in order to gr.nerate the parametric 
estimate. 
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22. A method according to claim 19, wherein one of the 
methods on the basis of which the parametric estimate is 
obtained ia chosen from the group consisting of: causal 
Simultaneous Auto- Regressive Moving Average (SARMA) models, 
non- causal SARMA models, Gaussian Markov Random Field models 
and Maximum Entropy Spectral Estimates and wherein the order 
of the model is identified, the parameters of the model are 
estimated and the spectrum of the estimated model is gener- 
ated. 

23. A method according to claim 22, wherein the method on the 
basis of which the parametric estimate is obtained, is a non- 
causal SARMA model using an estimator according to claim S8. 

24. A method according to any of the preceding claims, where- 
in at least one feature is related to parameters of a contour 
encompassing at least a predetermined percentage of the 
energy of the power spectrum. 

25. A method according to claim 24, wherein the contour is 
determined by 

(a) defining a contour, in terms of one axes in each 

dimension, around the centre of the power spectrum, 
all points on the contour having substantially the 
same distance to the centre of the power spectrum, 
and the contour encompassing less than the predeter- 
mined percentage of the energy of the power spectrum, 

<b) for each dimension of the data calculating the per- 
centage of the energy of the power spectrum 
encompassed by a dilated contour in which the axis in 
the dimension in question is increased by a predeter- 
mined distance, 

(c) increasing the axis of the contour, in the dimension 
in which the increase in energy encompassed by the 
dilated contour is the largest, by the predetermined 
distance and 
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(d) repeating steps (b) and (c) until the percentage 
encompassed by the contour exceeds or equals the 
predetermined percentage . 

26. A method according to claim 24 or 25, wherein the contour 
5 is ellipsoidal, optionally spherical. 

27. A method according to claim 25 or 26, wherein the axis of 
the contour are orthogonal and are defined along principal 
directions of the power spectrum. 

28. A method according to claim 1, wherein procedure (c) is 
10 used and wherein the Fourier method comprises 

subjecting the image data to Fourier transformation, 
performing a subsequent mathematical transformation of 
the Fourier transformed data, 

converting the transformed data back into the spatial 
15 domain. 

29. A method according to claim 28, wherein the subsequent 
mathematical transformation is linear or non- linear. 

30. A method according to claim 29, wherein the subsequent 
mathematical transformation is performed by raising the data 

20 to a power larger than l, such as a power in the range of 
1.1-10, preferably 1.5-4, such as about 2.0. 

31. A method according to any of claims 28-30, wherein a 
third mathematical transformation is performed in order to 
preserve only part of the magnitude information of the trans - 

25 formed data. 

32. A method according to any of claims 28-31, wherein the 
emphasized information in the image data relating to the 
trabecular structure is subjected to a grey- scale morphologi- 
cal operation. 
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33. A method according to claim 32, wherein the grey-scale 
morphological operation is a top- hat operation wherein all 
pixels fitting in the top-hat are assigned a value different 
from a specified background pixel value. 

5 34. A method according to claim 33, wherein the image data 
which has been subjected to the top- hat operation is thres- 
holded so as to assign a foreground pixel value to the pixels 
having a pixel value different from the background pixel 
value . 

10 35. A method according to any of claims 32-34, wherein the 
image data which has been subjected to the grey- scale mor- 
phological operation is subjected to further noise reduction 
to substantially remove isolated single pixels and pixel 
groups smaller than a predetermined threshold number of 

15 pixels having a pixel value different from the background 
pixel value. 

36. A method according to any of claims 32-35, wherein an 
operation is performed comprising determining, for each 
background pixel, the distance according to a given metric 

20 from the background pixel to the nearest foreground pixel, 

and wherein at least one feature is extracted from the result 
of the operation. 

37. A method according to claim 36, wherein at least one 
feature relates to a mean value and/or standard deviation 

25 and/or the coefficient of variation and/or skewness and/or 
kurtosis of the determined distances from all background 
pixels to the nearest foreground pixel. 

38. A method according to any of claims 32-35, wherein an 
operation is performed comprising determining a measure for 

30 each background pixel based upon determining the distance 

from the background pixel in a number of given directions in 
the image data to the nearest foreground pixel and at least 
one feature is extracted from the result of the operation. 
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39. A method according to claim 38, wherein at least one 
feature relates to a mean value and/or standard deviation 
and/or the coefficient of variation and/or skewness and/or 
kurtosis of the determined measures for all background 

5 pixels . 

40. A method according to claim 38 or 39, wherein the measure 
for each background pixel is a mean value and/or standard 
deviation and/or the coefficient of variation and/or skewness 
and/or kurtosis and/or maximum distance of the determined 

10 distances in the given directions. 

41. A method according to any of the preceding claims, 
wherein at least one feature related to the bone or to the 
vertebrate, such as age and/or sex, and/or species, and/or 
race and/or the specific bone considered in the vertebrate, 

15 and/or an estimated Bone Mineral Density of the bone, and/or 
an estimated Bone Mineral Content of the bone, is included in 
the estimation procedure. 

42. A method according to claim 41, wherein the Bone Mineral 
Density is estimated by including data from a reference 

20 object in the exposure of the bone to the electromagnetic 

radiation and where the Bone Mineral Density is estimated by 
correlating the absorption of the electromagnetic radiation 
of the bone and of the reference object. 

43. A method according to any of the preceding claims, 

25 wherein the estimation procedure is based on a statistical 
model, taking into account the correlation structure in the 
data set, so as to assign appropriate weights to the 
significant features in accordance with the predetermined 
relationship. 

30 44. A method according to any of the preceding claims, 

wherein the reference bone quality is an absolute or relative 
bone quality of bones for which image data have been: 
processed in accordance with any of the preceding claims. 
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45. A method according to claim 44, wherein the reference 
bone quality is determined on the basis of measurement of the 
mechanical bone strength, and/or Bone Mineral Density 
measurement, and/or on the basis of a Bone Mineral Content 

5 measurement, and/or on the basis of a score of the bone 
strength by a skilled radiologist. 

46. A method according to any of the preceding claims, 
wherein the predetermined relationship between the features 
and reference bone quality parameters is established on the 

10 basis of reference bone quality parameters and features 
extracted according to any of the preceding claims. 

47. A method according to claim 46, wherein the predetermined 
relationship is defined in terms of a model taken from the 
group consisting of: a General Linear Model, a Generalized 

15 Linear Model, an Artificial Neural Network, a Causal 

Probabilistic Net and Classification And Regression Trees. 

48. A method according to any of the preceding claims, 
wherein the vertebrate is a human. 

49. A method according to any of the preceding claims, 
20 wherein the vertebrate is a horse. 

50. A method according to any of the preceding claims, 
wherein the vertebrate is a large ape, a great ape or an 
anthropoid ape. 

51. A method according to any of the preceding claims, 
25 wherein the vertebrate is a pig. 

52. A method according to any of the preceding claims, 
wherein the vertebrate is a cow. 

53. A method according to any of the preceding claims, 
wherein the bone is taken from the group consisting of 

30 radius, femur, corpus vertebrae (LI, L2, L3, L4, L5, Tl, T2, 
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T3, T4, T5, T6, T7, T8, T9 , T10, Til, T12, CI, C2, C3, C4, 
C5, C6 f C7) , calcaneus, talus, os carpi, metatars, metacarpi, 
falanges, tibia, fibula, patella, ulna, humerus, mandible, 
clavicula, scapula, os coxae, os naviculare, os cuboideum, os 
5 cuneiform I, os cuneiform II, os cuneiform III, os sacrum, os 
coccygis . 

54. A method according to claim l, wherein procedure (d) is 
used and wherein the features relating to the periodicity of 
the trabecular structure of the part of the bone are derived 
10 from information relating to the height, width and/or 

curvature of peaks of and/or to the total area under peaks of 
the power spectrum of the image data, the smoothed 
periodogram of the image data, or the estimate of the 
parametric estimate of the power spectrum of the image data. 

15 55. A method according to claim 54, wherein the information 
is derived relating to the height, width and/or curvature of 
peaks of the smoothed periodogram of the image data and/or 
total area under the peaks thereof, and wherein the position 
of the peaks are identified in the power spectrum of the 

20 image data. 

56. A method for estimating the bone quality of a vertebrate, 
on the basis of two-dimensional image data comprising 
information relating to the trabecular structure of at least 
a part of a bone of the vertebrate, the image data being data 
25 obtained by exposing at least the part of the bone to 

electromagnetic radiation, the method comprising subjecting 
the image data to a statistical analysis comprising 



a background correction procedure in which low frequent 
intensity variations not related to the trabecular 
30 structure of the bone is reduced relative to image data 

related to the trabecular structure of the part of the 
bone, 
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an image manipulation and feature extraction procedure 
comprising subjecting the image data to at least one of 
the following procedures: 

(a) obtaining an estimate of the parametric estimate 
5 of the power spectrum of the image data and 

extracting features relating to the energy 
distribution of the parametric estimate, 

(b) obtaining, on the basis of image data on which a 
Fourier method has been used to emphasize the 

10 information in the image data relating to the 

trabecular structure, an estimate of a grey- level 
co-occurrence matrix and extracting at least one 
feature on the basis of the estimated co- 
occurrence matrix, 
15 (c) obtaining an estimate of the projected trabecular 

pattern of the image data by using a Fourier 
method to emphasize the information in the image 
data relating to the trabecular structure and 
subjecting the manipulated image data to a 
20 morphological operation, and extracting features 

relating to the trabecular structure from the 
estimated projected trabecular pattern, 
(d) obtaining, on the basis of a frequency analysis 
the power spectrum of the image data, a smoothed 
25 periodogram of the image data, or a parametric 

estimate of the power spectrum of the image data, 
features relating to the periodicity of the 
trabecular structure of the part of the bone, 
and an estimation procedure in which the bone quality of 
30 the vertebrate is estimated with a Multiple Correlation 

Coefficient better than 0.5 on the basis of the derived 
features and optionally other features related to the 
bone or the vertebrate and a predetermined relationship 
between the features and reference bone quality 
35 parameters. 

S7. A method according to claim 56, wherein the bone quality 
of the vertebrate is estimated with a Multiple Correlation 
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Coefficient better than 0.55, such as better than 0.6, 
preferably better than 0.65, such as better than 0.7, 
preferably better than 0.8, such as better than 0.85, 

58. A method for estimating Non-Causal Simultaneous Auto- 

5 Regressive Moving Average models in two or more dimensions, 
the method comprising optimizing a given direct measure of 
the flatness of the residual spectrum of the model. 

59. A method according to claim 50, wherein the optimization 
of the given flatness measure is obtained by 

10 (a) generating a set of initial parameters for the model, 
(b) generating the residual spectrum of the model on the 



2 0 60. A method according to claims 58 or 5, wherein the measure 
of flatness of the residual spectrum is obtained as the 
entropy of the normalized residual power spectrum, dis- 
regarding the DC value. 



15 



(c) 



<e) 



(d) 



basis of the parameters, 

obtaining the measure of the flatness of the residual 
spectrum, 

obtaining a new iterate of the parameters on the 
basis of the flatness measure and a search direction 
in parameter space, 

repeating steps {b) - (d) until given stop criterion is 
reached. 
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Fig. 6 
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Fig. 10 
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Fig. 12 
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Exact torus-ML estimates vs. true parameters (isotropic case) 
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MORSE estimates vs. true parameters (isotropic case) 
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